Computational aspects of astrophysical MHD and turbulence 



o 



Axel Brandenburg 
Nordita, Blegdamsvej 17, DK 2100 Copenhagen 0, Denmark; 
Department of Mathematics, University of Newcastle upon Tyne NEl 7RU, UK 

o : 

Abstract 

The advantages of high-order finite difi'erence scheme for astrophysical MHD and turbulence sim- 
ulations are highlighted. A number of one-dimensional test cases are presented ranging from various 

■ shock tests to Parker-type wind solutions. Applications to magnetized accretion discs and their as- 
^0 ' sociated outflows are discussed. Particular emphasis is placed on the possibility of dynamo action in 

, three-dimensional turbulent convection and shear flows, which is relevant to stars and astrophysical 

■ discs. The generation of large scale flelds is discussed in terms of an inverse magnetic cascade and 
the consequences imposed by magnetic helicity conservation are reviewed with particular emphasis 

, on the issue of a-quenching. 

> ■ 

0\ '. 1 Introduction 
^ ■ 

Over the past 20 years multidimensional astrophysical gas simulations have become a primary tool to 
understand the formation, evolution, and the final fate of stars, galaxies, and their surrounding medium. 
The assumption that those processes happen smoothly and in a non-turbulent manner can at best be 
regarded as a first approximation. This is evidenced by the ever improving quality of direct imaging 
techniques using space telescopes for example. At the same time not only have computers become large 
enough to run three-dimensional simulations with relatively little effort, there have also been substan- 
O I tial improvements in the algorithms that are used. In fact, there is now a vast literature on numerical 

-4— > ' astrophysics. An excellent book was published recently by LeVeque et al. (1998) where both numerical 

^ \ methods and astrophysical applications were discussed in great detail. Most of the applications focused 

■ however on rather more 'violent' processes such as supersonic jets, supernova explosions, core collapse, 
. 5^ \ and on radiative transfer problems, while hydromagnetic phenomena and turbulence problems where only 

■ touched upon briefly. Meanwhile, hydromagnetic turbulence simulations have become crucial for under- 
\ standing viscous dissipation in accretion discs (Hawley, Gammie, & Balbus 1995), and for understanding 

■ - - I magnetic field generation by dynamo action in discs (Brandenburg et al. 1995, 1996a, Hawley, Gammie, 

& Balbus 1996, Stone et al. 1996), stars (Nordlund et al. 1992, Brandenburg et al. 1996b), and planets 
(Glatzmaier & Roberts 1995, 1996). 

Much of the present day astrophysical hydrodynamic work is based on the ZEUS code, which has 
been documented in great detail and described with a number of test cases in a series of papers by 
Stone & Norman ( 1992a, b). The main advantage is its flexibility in dealing with arbitrary orthogonal 
coordinates which makes the code applicable to a wide variety of astrophysical systems. The code, which 
is freely available on the net, uses artificial viscosity for stability and shock capturing, and is based on 
an operator split method with second-order finite differences on a staggered mesh. Another approach 
used predominantly in turbulence research are spectral methods (e.g. Canuto et al. 1988), which have 
the advantage of possessing high accuracy. Although these methods are most suitable for incompressible 
flows (imposing the solenoidality condition is then straightforward), they have also been applied to 
compressible flows (e.g. Passot & Pouquet 1987). As a compromise one may resort to high order finite 
difference methods, which have the advantage of being easy to implement and yet have high accuracy. 
Compact methods (e.g. Lele 1992) are a special variety of high order finite difference methods, but the 
truncation error is smaller than for an explicit scheme of the same order. Compact schemes have been 
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used by Nordlund & Stein (1990) in simulations of solar convection (Stein & Nordlund 1989, 1998) and 
convectivc dynamos (Nordlund et al. 1992, Brandenburg et al. 1996b), for example. 

The use of compact methods involves solving tridiagonal matrix equations, making this method es- 
sentially nonlocal in that all points are now coupled at once. This is problematic for massively parallel 
computations, which is why Nordlund & Galsgaard (1995, see also Nordlund, Galsgaard, & Stein 1994) 
began to use explicit high order schemes for their work on coronal heating by reconnection (Galsgaard 
& Nordlund 1996, 1997a, b). In their code the equations are solved in a semi-conservative fashion using 
a staggered mesh. This code was also used by Padoan, Nordlund, & Jones (1997) and Padoan, & Nord- 
lund (1999) in models of isothermal interstellar turbulence in molecular clouds, and by Rognvaldsson, 
Nordlund, & Sommer-Larsen (2001) in simulations of cooling flows and galaxy formation. 

A somewhat different code was used by Brandenburg (1999) and Bigazzi (1999) in simulations of 
the inverse magnetic cascade, by Kerr & Brandenburg (1999) in work on the possibility of a singularity 
of the nonresistive and inviscid MHD equations, and by Sanchez-Salcedo & Brandenburg (1999, 2001) 
in simulations of dynamical friction. A two-dimensional version of the code modelling outflows from 
magnetized accretion discs has been described by Brandenburg et al. (2000). This code uses sixth order 
explicit finite differences in space and third order Runge-Kutta timestepping. It employs central finite 
differences, so the extra cost of recentering a large number of variables between staggered meshes each 
timcstcp is avoided. 

Apart from high numerical accuracy, another important requirement for astrophysical gas simulations 
is the capability to deal with a large dynamical range in density and temperature. This requirement 
favors the use of non-conservative schemes, because then logarithmic variables can be used which vary 
much less than linear density and energy density per unit volume. Solving the nonconservative form of 
the equations can be more accurate than solving the conservative form. The conservation properties can 
then be used as an indicator for the overall accuracy. 

In this chapter we concentrate on numerical astrophysical turbulence aspects starting with a discussion 
of different numerical methods and a description of the results of various numerical test problems. This 
is a good way of assessing the quality of a numerical scheme and of comparing with other methods; see 
Stone & Norman (1992a,b) for a series of tests using the ZEUS code. After that we discuss particular 
astrophysical applications including stellar convection, accretion disc turbulence and associated outflows, 
as well as the generation of magnetic fields (small scale and large scale) from turbulence in various 
astrophysical settings. 



2 The Navier-Stokes equations 

The discussion of magnetic fields will be postponed until later, because the inclusion of the Lorentz force 
in the momentum equation is straightforward. We begin by writing down the Navier-Stokes equations in 
nonconservative form and rewrite them such that the main thermodynamical variables are entropy and 
either logarithmic density or potential enthalpy. These variables have the advantage of varying spatially 
much less than for example linear pressure and density. 
The primitive form of the continuity equation is 

| = -v.(p.), (1) 

which means that the local change of density is given by the divergence of the mass flux at that point. 
The Navier-Stokes equation can be written as 

p—=-Vp-pV<^ + F + V-T, (2) 

where D/Dt = d / dt + u ■ 'V is the advective derivative, p is the pressure, $ is the gravitational potential, 
F is a body force (e.g., the Lorentz force), and r is the stress tensor. 

The Navier-Stokes equation is here written in terms of forces per unit volume. As argued above, if 
the density contrast is large it is advantageous to write it in terms of forces per unit mass and to divide 
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by p. Before we can replace p and p by entropy and logarithmic density or potential entropy we first have 
to define some thermodynamic quantities. 

Internal energy, e, and specific enthalpy, h are related to each other by 

h = e+ pv, (3) 

where v = 1/p is the specific volume and p the density. The specific entropy is defined by 

Tds^de+pdv, (4) 

where T is temperature. The specific heats at constant pressure and constant volume are defined as 
Cp = dh/dT\p and c„ = de/dT|„, their ratio is 7 = Cp/cy, and their difference is TZ/p = Cp — c„, where TZ 
is the universal gas constant and p the specific molecular weight. 

In the following we assume Cp and c„ to be constant for all processes considered. Ionization and 
recombination processes are therefore ignored here, although this is not a major obstacle; see, e.g., 
simulations of Nordlund (1982, 1985), Steffen, Ludwig, & Kriifi (1989), Stein & Nordlund (1989, 1998), 
Rast et al. (1993), and Rast & Toomre (1993a, b) where realistic equations of state have been used. 

We now assume that Cp and Cy are constant, so internal energy and specific enthalpy are given by 

h — CpT and e — c^T, (5) 
This allows us then to write the specific entropy (up to an additive constant) as 

s — Cylnp ~ Cp In p. (6) 
The pressure gradient term in the momentum equation can then be written as 

-Vp ^ -V Inp ^ ^{V In p + Vs/cp) ^ cliV In p + Vs/cp), (7) 
P P P 

where we have used 

Cs = 7P/P = — exp[(7- l)ln(p/po) +7s/cp] = c^o f— ) exp(7s/cp), (8) 
Po \Po/ 

where Cg is the adiabatic sound speed, and c^q — jpo/po- 

With these preparations the evolution of velocity m, logarithmic density In p, and specific entropy s 
can be expressed as follows: 



^ = -c2(V Inp + Vs/cp) - V$ + / + i V • (2^pS), (9) 
Dlnp 



Bt 



= -V-M, (10) 



r^ = 2j.s2 + r-pA, (11) 



where f = F / p is the body force per unit mass, F and A are heating and cooling functions, v kinematic 
viscosity and S is the (traceless) strain tensor with the components 

Stj = |(wij + - lSi-jUk,k)- (12) 

In the presence of an additional kinematic bulk viscosity, ^, the term 2iySij under the divergence in eq. (^ 
would need to be replaced by 2i/Sij + C^ij^ ' and the viscous heating term, 2vS'^, in eq. ( pl| ) would 
need to be replaced by 2i/S^ + C(V • m)^. 
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Instead of using p as a dependent variable an can also use the specific enthalpy ft., which allows us to 
write the pressure gradient as 

--Vp= -Vh + TVs. (13) 
P 

This formulation is particularly useful if the entropy is nearly constant (or if the gas is barotropic, i.e. 
p — p{p)) and if there is a gravitational potential so that the potential enthalpy H = h + $ can be 
used as dependent variable. In order to express eq. (|l^) in terms of h we write down the total differential 
of the specific entropy, 

ds/c„ = -dlnp - dlnp = -dln/i - ( 1 - ^ ) dlnp, (14) 
7 7 V 7/ 

so 

Dln/i Ds/c„ , ,Dlnp , , 
=7 — ^ + (7-1 -■ 15 

Furthermore, TVs = h'^s/cp, and so the final set of equations is 









— = -VH + hV.s/cp + f- 


h V • (2j.pS), 


(16) 


Dt h ^ 


-pA), 


(17) 


„^ , Ds/cp 
Dt =-V<f + 7^ Dt 




(18) 


- V • u, 



where we have absorbed $ in the potential enthalpy H — h + ^. In this formulation the density can be 
recovered as 



Po 



1/(7-1) 

(19) 



7/ \cto, 

(in dimensional form) or, for 7 = 5/3 and in nondimensional form (where pg = pg = Cp = 1), 

p= (0.4/^)l■5e~2■^^ (20) 

We shall use either of the two sets of the equations, (^-(pl]) or (|l^)-(|l^), in some of the following 
sections, especially in connection with shock tests and stellar wind problems. In these cases the gravity 
potential $ is important and it turns out that the potential enthalpy H = h + ^ varies only very little 
near the central object even though $ itself tends to become singular. 

The heating and cooling terms (F and A) are important for example in the case of interstellar tur- 
bulence which is driven primarily by supernova explosions which inject a certain amount of thermal 
energy (J pTdV) with each supernova explosion. MHD turbulence simulations of this type were per- 
formed recently by Korpi et al. (1999). At the same time there is cooling through various processes 
(e.g. bremsstrahlung at high temperatures) which transports energy either nonlocally via a cooling term 
A(T), or locally via thermal conduction or radiative diffusion. In the radiative diffusion approximation 
we express A as —p^A = V • KWT , where K is the radiative conductivity which is in general a function 
of temperature and density. The radiative diffusivity (which has the same dimensions as the kinematic 
viscosity v) is given by x = K/{pCp)^ so 

- pk/{cpT) = V . pxCpVT. (21) 

Since we shall use a nonconservative scheme with centered finite differences is it important to isolate 
second derivative terms, so 

- pA//i = x(V^lnT + Vlnp- VlnT), (22) 
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where we have assumed for simplicity that x is constant. In terms of s/cp and hip we have 

- pA/h = X7 [V^s/cp + Vad In p + 7( Vs/cp + V In p) • ( Vs/cp + Vad V \n p)] , (23) 

where Vad = 1 — l/7isa commonly used abbreviation in stellar astrophysics. For 7 = 5/3 we have 
Vad = 2/5 = 0.4. We shall use eq. ( |2^ ) later in connection with shock and wind calculations. However, 
we begin by discussing first a suitable numerical scheme which will be used in most of the cases presented 
below. 

3 The advantage of higher-order derivative schemes 

Spectral methods are commonly used in almost all studies of ordinary (usually incompressible) turbu- 
lence. The use of this method is justified mainly by the high numerical accuracy of spectral schemes. 
Alternatively, one may use high order finite differences that are faster to compute and that can possess 
almost spectral accuracy. Nordlund & Stein (1990) and Brandenburg et al. (1995) use high order finite 
difference methods, for example fourth and sixth order compact schemes (Lele 1992). [| 

In this section we demonstrate, using simple test problems, some of the advantages of high order 
schemes. We begin by defining various schemes including their truncation errors and their high wavenum- 
ber characteristics. We consider centered finite differences of 2nd, 4th, 6th, 8th, and 10th order, which 



are given respectively by the formulae 

/; = (-/.-i+/.+i)/(25x), (24) 

n - {h-2 - 8/,-i + 8/,+i - /,+2)/(12fa), (25) 

/; = + 9/^-2 - 45/,;_i + 45/,+i - 9/,+2 + /,+3)/(60fe), (26) 

/; = (3/,_4 - 32/,_3 + 168/,_2 - 672/,_i + 672/,+i - 168/,+2 + 32/,;+3 - 3/^+4) /(8405a:), (27) 

= (-2/,_5 + 25/,_4 - 150/,_3 + 600/,_2 - 2100/,_i + 2100/,+i - ...)/(25205a;), (28) 
for the first derivative, and 

= (/,_i - 2u + .n+i)/{Sx'), (29) 

f- = {-h-2 + m-i - 30/, + 16/,+i - /,+2)/(12fe2), (30) 

= (2/,_3 - 27/,_2 + 270/,_i - 490/, + 270/,+i - 27/,+2 + 2/,+3)/(180fe2), (31) 



/f = (-9/,_4 + 128/,;_3 - 1008/,_2 + 8064/,_i - 14350/, + 8064/,+i - 1008/,+2 + ...)/(5040fe2), (32) 

f- - (8/,_5 - 125/,_4 + 1000/,;_3 - 6000/,_2 + 42000/,;_i - 73766/, + 42000/,+i - ...)/i252006x^), (33) 

for the second derivative. The expressions for one-sided and semi-onesided finite difference formulae are 
given in Appendix 

3.1 High wavenumber characteristics 

The chief advantage of high order schemes is their high fidelity at high wavenumber. Suppose we differ- 
entiate the function sin /ex, we are supposed to get kcoskx, but when k is close to the Nyquist frequency, 
^Ny = n/Sx, where Sx is the mesh spacing, numerical schemes yield effective wavenumbers, k^g, that can 
be significantly less than the actual wavenumber k. Here we calculate k^s from 

(cos kx)[^^^ = -fccff sin kx. (34) 

^The fourth order compact scheme is really identical to calculating derivatives from a cubic spline, as was done in 
Nordlund & Stein (1990). In the book by CoUatz (1966) the compact methods are also referred to as Hermitian methods 
or as Mehrstellen-Verfahren, because the derivative in one point is calculated using the derivatives in neighboring points. 
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When k = fcNy, every centered difference scfieme will give kcs = 0, because then the function values 
of cos fea; are just — 1,+1,— 1, so the function values on the left and the right are the same, and the 
difference that enters the scheme gives therefore zero. 

If is useful to mention at this point that for a staggered mesh, where the first derivative is evaluated 
between mesh points, the value of the first derivative remains finite at the Nyquist frequency, provided 
one does not need to remesh back to the original mesh. Especially in the context with magnetic fields, 
however, remeshing needs to be done quite frequently, which therefore diminishes the advantage of a 
staggered mesh. 

In figure ^ we plot effective wavenumbers for different schemes. Apart from the different explicit finite 
difference schemes given above, we also consider a compact scheme of 6th order, which can be written in 
the form 

+ fl + = (/-2 - 28/,_i + 28/,+i - /,+2)/(36fe), (35) 

for the first derivative, and 

Tifr'-i + /" + TTf"+i = (3/^-2 + m-i - 102/, + 48/,+i + 3,U+2)/{USx^). (36) 

for the second derivative. As we have already mentioned in the introduction, this scheme involves obvi- 
ously solving tridiagonal matrix equations and is therefore effectively nonlocal. 



1st derivative 2nd derivative 




Figure 1: Effective wave numbers for first and second derivatives using different schemes. Note that for 
the second derivatives the sixth order compact scheme is almost equivalent to the tenth order explicit 
scheme. For the first derivative the sixth order compact scheme is still superior to the tenth order explicit 
scheme. 

In the second panel of figure |l| we have plotted effective wavenumbers for second derivatives, which 
were calculated as 

(cos fca;)(j„,jj = — fcoff cos kx. (37) 

Of particular interest is the behavior of the second derivative at the Nyquist frequency, because that 
is relevant for damping zig-zag modes. For a second-order finite difference scheme k'^g is only 4, which 
is less than half the theoretical value of tt^ = 9.87. For fourth, sixth, and tenth order schemes this 
value is respectively 5.33, 6.04, 6.83. The last value is almost the same as for the 6th order compact 
scheme, which is 6.86. Significantly stronger damping at the Nyquist frequency can be obtained by using 
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hyperviscosity, which Nordlund & Galsgaard (1995) treat as a quenching factor that diminishes the value 
of the second derivative for wavenumbers that are smaU compared with the Nyquist frequency. Accurate 
high order second derivatives (with no quenching factors) are important when calculating the current 
J in the Lorentz force J x B from a vector potential A using —fiaJ = W^A — W • A. This will be 
important in the MHD calculations presented below. 



3.2 The truncation error 

One can express /i-i, /i+i, etc, in terms of the derivatives of / at point i, so 

- /, - fe/; + ^Sx'll' - hx^ir + ... (38) 

= /. + sxn + Isx'f- + ifavr + - (39) 

Inserting this into the finite difference expressions yields for the second order formula 

(/02„d ^ - h-i)/{2Sx) = ,n + Isx^n"- (40) 

The error scales quadratically with the mesh size, which is why the method is called second order. The 
truncation error is proportional to the third derivative of the function. Because this is an odd derivative 
it corresponds to a dispersive (as opposed to diffusive) error. Schemes that are only first order (or of any 
odd order) have diffusive errors, and it is this what is sometimes referred to as numerical diffusivity, which 
is not to be confused with artificial diffusivity that is sometimes used for stability and shock capturing. 
For the other schemes given in eqs. (|25|)-(p^) the truncation errors are 

(/;)4th = /;+3x io-2faVi'\ (41) 

(/06th = /; + 7x lO-^feVf", (42) 

(/;)ioth = /; + 3x 10-4<5a;i"/f^ (43) 

For the sixth order compact scheme the error scales like for the sixth order explicit scheme, but the 
coefRcient in front of the truncation error is about ten times smaller, so 

(^,)Con.pact ^ ^, ^ ^ 1Q-Ux'ft'''''>. (44) 

For the second derivatives we have 

(/n2nd = /r+8xio-2fevf\ (45) 

(/n4th = /r-ixio-2favf'\ (46) 
(/r)6th = /;'+2xio-3fevr\ (47) 

(/nioth-/r-5xl0-5fei«/f). (48) 

Again, for the sixth order compact scheme the scaling is the same as for the sixth order explicit scheme, 
but the coefficient in front of the truncation error is about 5 times less, so 

(/r)6~* = /r + 3.2 X 10-4 Sx^ft'''\ (49) 

This information about the accuracy of schemes would obviously be of little use if the various schemes 
did not perform well when applied to real problems. For this reason we now begin by carrying out various 
tests, including advection and shock tests. 
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3.3 Advection tests 



As a first test we compare the various schemes by performing inviscid advection tests and solve the 
equation Df /Dt — 0, i.e. 

/ = -uf\ (50) 

on a periodic mesh. It is advantageous to use a relatively small number of meshpoints (here we use 
Nx = 8 meshpoints), because that way we see deficiencies most clearly. This case is actually also relevant 
to real applications, because in practice one will always have small scale structures that are just barely 
resolved. 

After some time an initially sinusoidal signal will suffer a change in amplitude and phase. We have 
calculated the amplitude and phase errors for schemes of different spatial order. For the time integration 
we use high order Runge-Kutta methods of 3rd or 4th order, RK3 and RK4, respectively. In most cases 
considered below we use the RK3 scheme that allows reasonable use of storage. It can be written in three 
steps (Rogallo 1981) 

1st step: / = / + liStf, g = f + C,5tf, 

2nd step: f = g + j2Stf, g^f + C2Stf, (51) 
3rd step: f ^ 9 + Iz^tf, 

where 

8 5 3 17 5 

^^=15' ^^^12' ^^=4' ^^ = -60' ^' = --2- 

where, / and g always refer to the current value (so the same space in memory can be used), but / is 
evaluated only once at the beginning of each of the three steps at < t^, ti/^ = to + li5t w to + 0.5333(5t, 
and at t2/3 — t-o + {li + Ci + l2)5t = + '^St. Even more memory-effective are the so-called 2iV-schemes 
that require one set of variables less to be hold in memory. Such schemes work for arbitrarily high order, 
although not all Runge-Kutta schemes can be written as 2iV-schemes (Williamson 1980, Stanescu & 
Habashi 1998). These schemes work iteratively according to the formula 

Wi = aiWi-i+ 5tF{U^i,Ui^i), Ui ^ Ui^i + PiWi. (53) 

For a three-step scheme we have « = 1, ...,3. In order to advance the variable u from u'-"-' at time t^"-* to 
at time = + 5h we set in eq. (||) 



= and = U3, (54) 



with ui and U2 being intermediate steps. In order to be able to calculate the first step, i = 1, for which 
no Wi-i = Wo exists, we have to require ai = 0. Thus, we are left with 5 unknowns, a2, as, (32, and 
P3. Three conditions follow from the fact that the scheme be third order, so we have to have two more 
conditions. One possibility is the choose the fractional times at which the right hand side is evaluated, for 
example (0, 1/3, 2/3) or even (0, 1/2, 1). In the latter case the right hand side is evaluated twice at the 
same time. It is therefore some sort of 'predictor-corrector' scheme. In the following these two schemes 
are therefore referred to as 'symmetric' and 'predictor/corrector' schemes. Yet another possibility is to 
require that inhomogeneous equations of the form -it = t" with n = 1 and 2 are solved exactly. Such 
schemes are abbreviated as 'inhomogeneous' schemes. The detailed method of calculating the coefficients 
for such third order Runge-Kutta schemes with 27V-storage is discussed in detail in Appendix |^. Several 
possible sets of coefficients are listed in table |^ and compared with the favorite scheme of Williamson 
(1980). Note that the first order Euler scheme corresponds to /3i = 1 and the classic second order to 
a2 = -1/2, I3i = 1/2, and /32 = 1. 

We estimate the accuracy of these schemes by solving the homogeneous differential equation 



-l/n 



m(1) = 1. (55) 



The exact solution is m = V\ In table ^ we list the rms error with respect to the exact solution, for the 
range 1 < t < 4 and fixed time step 6t = 0.1 using n = —1, 2 or 3. 
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Table 1: Possible coefficients for different 2iV-RK3 schemes. 



label 


a2 


"3 


Pi 


P2 


/33 


(ti - tn)/6t 


[h - h)l5t 


symmetric (i) 


-2/3 


-1 


1/3 


1 


1/2 


1/3 


2/3 


symmetric (ii) 


-1/3 


-1 


1/3 


1/2 


1 


1/3 


2/3 


predictor /corrector 


-1/4 


-4/3 


1/2 


2/3 


1/2 


1/2 


1 


inliomogeneous 


-17/32 


-32/27 


1/4 


8/9 


3/4 


1/4 


2/3 


quadratic 


-0.367 


-1.028 


0.308 


0.540 


1 


0.308 


0.650 


Williamson (1980) 


-5/9 


-153/128 


1/3 


15/16 


8/15 


1/3 


3/4 



Table 2: Errors (in units of 10 ^) for different 2iV-RK3 schemes, obtained by solving eq. ( ^5| ) in the 
range 1 <t < A with 5t = 0.1 and different values of n. 



label 


n = -1 


n = 2 


n ~ i 


symmetric (i) 


69 


103 


193 


symmetric (ii) 


226 


119 


411 


predictor/corrector 


469 


346 


1068 


inhomogeneous 


84 


6 


97 


quadratic 


197 


94 


339 


Wilhamson (1980) 


68 


10 


123 


for comparison: RK3 


66 


13 


134 



The length of the time step must always be a certain fraction of the Courant-Friedrich-Levy condition, 
i.e. 5t = kcFhSx/Utnn^, where fccFL = C(l) and t/max is the maximum transport speed in the system 
(taking into account advection, sound waves, viscous transport, etc). Too long a time step can not only 
lead to instability, but it also increases the error. 

In table ^ we give amplitude and phase errors for the various schemes. The most important conclusion 
to be drawn from this is the fact that low order spatial schemes result in large phase errors. In the case 
of a second order scheme the phase error is 36° after a single passage of a barely resolved wave through 
a periodic mesh. Higher order schemes have easily a hundred times smaller phase errors. The amplitude 
error, on the other hand, is virtually not affected by the spatial order of the scheme. The amplitude error 
is mainly affected both by the temporal order of the scheme and by the length of the timestep; see also 
table ^. Therefore, high order schemes with low dissipation and dispersion are particularly important in 
computational acoustics (Stanescu & Habashi 1998). However, in applications to turbulence a certain 
amount of viscosity is always necessary. This would decrease the amplitude of the wave further and would 
eventually be even more important. (This additional viscosity could be the real one, an explicit artificial, 
or an implicit numerical viscosity that would result from the discretisation error or the numerical scheme; 



see §|3J.) 

A common criticism of high order schemes is their tendency to produce Gibbs phenomena (ripples) 
near discontinuities. Consequently one needs a small amount of diffusion to damp out the modes near 
the Nyquist frequency. Thus, one needs to replace eq. ( ^o|) by the equation 

/ - -uf + lyf. (56) 

The question is now how much diffusion is necessary, and how this depends on the spatial order of the 
scheme. 

A perfect step function would produce large start-up errors; it is better to use a smoothed profile, for 
example one of the form 

, , /cosa;\ , , 

f{x) = tanh (-J-) , (57) 
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Table 3: Amplitude and phase errors for inviscid advection of the function / = cos kx with k = 2tt and 
N = 8 meshpoints after t = 20, corresponding to 20 revolutions in a periodic mesh. The amplitude error 
is counted positive when the amplitude decreases. A positive phase error means that the solution lags 
behind the theoretical one. 



St scheme 


2nd order 


4th order 


6th order 


10th order 


spectral 


RK4 


1% 


2% 


2% 


2% 


2% 


est = 0.6 


720° 


87° 


13° 


2.9° 


2.8° 


RK3 


10% 


14% 


14% 


15% 


15% 


est = 0.4 


716° 


83° 


8° 


-2.1° 


-2.3° 


RK3 


4% 


6.3% 


6.6% 


6.6% 


6.6% 


est = 0.3 


717° 


84° 


10° 


-0.6° 


-0.8° 



Table 4: Dependence of the amplitude and phase errors on the length of the timestep and the scheme 
used for the timestep. In all cases spectral x-derivatives are used. 



St scheme 


RK3 


RK3 


RK4 


RK4 


RK4 


est 


0.3 


0.4 


0.4 


0.6 


1.0 


amplitude error 


6.6% 


15% 


0.3% 


2% 


21% 


phase error 


-0.8° 


-2.3° 


0.6° 


2.8° 


18° 



where 6x is the mesh width. For a periodic mesh of length L one would obviously use / = f{kx), where 
k = 2tt/L. In that case the step width would be k5x. In the following we consider a periodic domain of 
size Lx — \ with ~ 60 meshpoints, so we use k — 27: / Lx — 2t:. 

In figure]^ we plot the result of advecting the periodic step-like function, f{kx), over 5 wavelengths, 
corresponding to a time T = L/u. The goal is to find the minimum diffusion coefficient v necessary to 
avoid wiggles in the solution. In the first two panels one sees that for a 6th order scheme the diffusion 
coefficient has to be approximately v — 0.01 u5x. For v = 0.005 u5x there are still wiggles. For a 10th 
order scheme one can still use v — 0.005 u5x without producing wiggles, while for a spectral scheme of 
nearly infinite order one can go down to 0.002 u5x without any problems. 

We may thus conclude that all these schemes need some diffusion, but that the diffusion coefficient 
can be much reduced when the spatial order of the scheme is high. In that sense it is therefore not true 
that high order schemes are particularly vulnerable to Gibbs phenomena, but rather the contrary! 

In figure ^ we compare the corresponding results of advection tests for second and fourth order 
schemes with the sixth order scheme. It is evident that a second order scheme requires a relatively high 
diffusion coefficient, typically around v — 0.05 wJx, but this leads to rather unacceptable distortions of 
the original profile. (It may be noted that, if one uses at the same time a 1st order temporal scheme, 
which has antidiffusive properties, and a time step which is not too short, then the antidiffusive error of 
the timestep scheme would partially compensate the actual diffusion and one could reduce the value of 
v, but this would be a matter of tuning and hence not generally useful for arbitrary profiles.) 

3.4 Burgers equation 

In the special case where the velocity itself is being advected, i.e. f = u, eq. ( p6| ) turns into the Burgers 
equation, 

u=-uu' + vu". (58) 
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6th order, i' = 0.0lu6x 6th order, i'=0.005u6x 




Figure 2: Resulting profile after advecting a step-like function 5 times through the periodic mesh. The 
dots on the solid line give the location of the function values at the computed meshpoints and the dotted 
line gives the original profile. For the panels on the right hand side the diffusion coefficient is too small 
and the profile shows noticeable wiggles. Sx — 1/60. 



In one dimension there is an analytic solution for a kink, 

M = C/(^l-tanh^-^^ , (59) 

where 5 — v/U is the shock thickness (e.g., Dodd et al. 1982). Note that the amplitude of the kink is twice 
its propagation speed. Expressed in terms of the Reynolds number, Re = UL/v, we have(5 = LRe-^ (We 
note in passing that the dissipative cutoff scale in ordinary turbulence is somewhat larger; 5 — LRe"'^^'*.) 
In order to have a stationary shock we use the initial condition 

M = -C/tanh(a;/2(5). (60) 

In figure ^ we present numerical solutions using the sixth order explicit scheme with different values of the 
mesh Reynolds number, 5x U /v, which was varied by changing the value of v. Here we used Nx = 100 
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Figure 3: Like figure |2| but for low-order schemes. Again, for the panels on the right hand side the 
diffusion coefficient is too small and the profile shows noticeable wiggles. For the 2nd-order scheme one 
needs a viscosity of i/ = 0.05 u5x to prevent wiggles, but then the resulting distortion of the original 
profile becomes rather unacceptable. 5x — 1/60. 



meshpoints in the range — 1 < x < 1. Note that the overall error, defined here as max(|/ — /cxacti), 
decreases with decreasing mesh width like 5x^. 

The test cases considered so far were not directly related to the Navier-Stokes equation, which permits 
sound waves that can pile up to form shocks, for example. This will be considered in the next section. 

3.5 Shock tube tests 

A popular test problem for compressible codes is the shock tube problem of Sod (1951). On the one hand, 
one can assess the sharpness of the various fronts. On the other hand, and perhaps most importantly, it 
allows one to test the conversion of kinetic energy to thermal energy via viscous heating. 

In the following we use the formulation of the compressible Navier-Stokes equations in terms of entropy 
and enthalpy, (|l^)-(|lj). We use units where po = pp = Cp = 1 and adopt the abbreviations A = Inp (not 
to be confused with the cooling function used in § y) . In one dimension (with z/ = const) these equations 
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Figure 4: Numerical solution of the Burgers equation using the sixth order explicit scheme with eq. ( |60[ ) 
as initial condition. In the left hand panel the lines give the exact solution ( |60| ) and the dots give the 
numerical solution for the corresponding value of the mesh Reynolds number: 6x U/v = 0.5 (solid line) 
1.0 (dotted line), and 2.0 (dashed line). In the right hand panel the scaling of the error with 5x is shown. 



reduce to 

u = -uu' ~ c2(A' + s') + v{v!' + A'w'), (61) 

s = -ws' + (7-l)£;w'2/c2 + Q^, (62) 

A = -mA' - u', (63) 

where dots and primes refer respectively to time and space derivatives, Qs describes the change of entropy 
due to radiative diffusion, and A = Inp is the logarithmic density. In eqs. (|6l|)-(|63|) we have used the 
abbreviation 

Cs = IvIp = 7 exp[(7 - 1) A + 7s] (64) 

for the adiabatic sound speed squared, and v = '^v \s the effective viscosity for compressive motions. 
This 4/3 factor comes from the fact that in 1-D 

S = diag (65) 

and therefore ^ 

-V • (2i/pS) = \v{u^,-j,^ + {}ci\p)^xUx.x\^ (66) 

so — fw^, or IvS^ — vv?^ ^. In the radiative diffusion approximation we have Qg ~ — Acond/(cpT), 
and so eq. (|2^) gives in one dimension 

Qs = Xl\s" + VadA" + 7(A' + s')(s' + VadA')] (67) 

In figure |^ we show the solution for an initial density and pressure jump of 1:10 and the the viscosity is 
now V — 0.6(5xcs. In this case a small amount of thermal diffusion (with Prandtl number xjv = 0.05) 
has been adopted to remove wiggles in the entropy. 

For stronger shocks velocity and entropy excess increase; see figures ^ and|^, where the initial pressure 
jumps are 1:100 and 1:1000, respectively, and the viscosities are chosen \,o\>qv — 1.6(5x Cs and v = 2A5x Cg. 
(For the stationary shock problem considered below we also find that the viscosity must increase with the 
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Figure 5: Standard shock tube test with an initial pressure jump of 1:10 and v = O.QSx Cg and xl^ — 0.05. 
The sohd line indicates the analytic solution (in the limit ^ 0) and the dots the numerical solution. 
Note the small entropy excess on the right of the initial entropy discontinuity. 

Mach number and, moreover, that the two should be proportional to each other.) In the cases shown in 
figures 1^ and ^ we were able to put x — without getting any wiggles in s. However, in the case of strong 
shocks (pressure ratio 1:1000) the discrepancy between numerical and analytical solutions becomes quite 
noticeable. 
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Figure 6: Same as figure ||, but for an initial pressure jump of 1:100 and v = \SS8xc^. 
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Figure 7: Same as figure ^ but for an initial pressure jump of 1:1000 and ly = 2.46x0^. 
In many practical applications shocks occur only in a small portion of space. One can therefore reduce 



14 



the viscosity outside shocks or, conversely, use a smah viscosity everywhere except in the locations of 
shock, i.e. where the flow is convergent (negative divergent). This leads to the concept of an artificial 
(Neumann-Richtmyer) shock viscosity. 



i^shock = Cshockfe^ {{-'^ ■ u) + )n.n. 



(68) 



where (...)n.n. indicates averaging over nearest neighbors and the subscript + means that only the positive 
part is taken. 

The last panel in figures ^ and ^ shows quite clearly how the entropy increases behind the shock. 
This entropy increase is just a consequence of the viscous heating term, vu''^ /h. Without this term the 
solution would obviously be wrong everywhere behind the shock, especially when the shock is strong. 

A somewhat simpler situation is encountered with standing shocks. In figure || we give an example of a 
numerically determined solution at Ma^lOO. The agreement in the jump for the numerically determined 
solution (dotted line and dots) and the theoretical solution (solid lines) is very good, although the position 
of the jump has moved away somewhat from the initial location (x = 0), but this is merely a consequence 
of having used more-or-less arbitrarily a tanh profile to smooth the initial jumps. After some initial 
adjustment phase the profiles do indeed remain stationary. Note also, however, that the entropy profile 
is slightly shifted relative to the profiles of u and A. 

It is interesting to note that when solving the Rankine-Hugionot jump conditions for shocks one is 
allowed to use the inviscid equations provided they are written in conservative form. Sometimes one finds 
in the literature the inviscid Navier-Stokes equations written in nonconservative form. This is not strictly 
correct, because without viscosity there would be no viscous heating and hence no entropy increase behind 
the shock. Moreover, it is quite common to consider a polytropic equation of state, p = Kp^ . Again, 
in this case the entropy is constant, and so energy conservation is violated. Nevertheless, given that 
polytropic equations of state are often considered in astrophysics we consider this case in more detail in 
the next subsection. 



3.6 Polytropic and isothermal shocks 

For polytropes with p — Kp^ , but F 7^ 7 in general, we can write 

TK 



-Vh + hVs = - 

so we can introduce a pseudo enthalpy h as 

~ FA' 



-Vp 

P 



F - 1 



P 



r-i 



r-i 



F- 1' 



This is consistent with a fixed entropy dependence, where s only depends on p like 



s{p) = -lnip/p-y) = -ln{Kp^-^) 

7 7 



-InK 

7 



F 

1 

7 



Inp, 



(69) 



(70) 



(71) 



which implies that in the polytropic case eq. (p^ ) is discarded. In the adiabatic case, F = 7, entropy is 
constant. In the isothermal case, p — c^p, we have F — > 1, so entropy is not constant, but it varies only 
in direct relation to — In p and not as a consequence of viscous heating behind the shock. 

In deriving the Rankine-Hugionot jump conditions one uses the conservation of mass, momentum, 
and energy in a comoving frame, where the following three quantities are constants of motion: 

7 P 



J = pu, I — pu + p, E 



(72) 



7- Ip 

The values of these three constants can be calculated when all three variables, u, p, and p, are known on 
one side of the shock. For polytropic equations of state, with p = Kp"^, the energy equation is no longer 
used, so there are only the following two conserved quantities, 

,2 



J = pu, I = pu + Kp"' 



(73) 
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Figure 8: Example of a very strong standing shock with Ma=100. Note the relative shift of the position 
where s increases relative to where A = In p increases. The viscosity is chosen to he u = Ma x 5x. 



The dependence of the velocity, density, pressure, and entropy jumps on the upstream Mach number is 
plotted in figure ^ for the case 7 = 5/3 and compared with the polytropic case using F = 7. 

Note that the pressure jump, P2/P1, is almost independent of the value of 7 and does also not 
significantly depend on the polytropic assumption. 
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Figure 9: The dependence of the velocity, density, pressure, and entropy jumps on the upstream Mach 
number for the case 7 = 5/3 (soUd hne) and comparison with the polytropic case using F = 7 (dotted 
hue). Note that the velocity and density jumps saturate at 1/4 and 4, respectively, while there is no such 
saturation for the polytropic shock. 




Figure 10: Like figure ^, but for 7 = F = 1.0001. Note that the velocity and density jumps saturate at 
much more extreme values than for 7 = 5/3. Thus, the results using polytropic and perfect gas equations 
agree up to much larger Mach numbers than for 7 = 5/3. 

4 Nonuniform and lagrangian meshes 

In many cases it is useful to consider nonuniform meshes, either by adding more points in places where 
large gradients are expected, or by letting the points move with the flow (lagrangian mesh). The la- 
grangian mesh is particularly useful in one-dimensional cases, because then the mesh topology (i.e. the 
ordering of mesh points) remains unchanged. Another method that gains constantly in popularity is 
adaptive mesh refinement (e.g., Grauer, Marliani, & Germaschewski 1998), which will not be discussed 
here. 

4.1 Nonuniform topologically cartesian meshes 

Nonuniform meshes can be implemented relative easily when each of the new coordinates depend on 
only one variable, for example when x — a;(x), y — y{y), and z — z{z). Here, x, y and z are cartesian 
coordinates on a uniform mesh. In the more general case, however, we have 

x^x{x,y,z), y = y{x,y,z), z^z{x,y,z), (74) 
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so that X, y and z derivatives of a function / can be calculated using the chain rule, 



df df dx 
dx dx dx 



df dy df dz _ df dxi 
dy dx dz dx dxi dx 



Corresponding formulae apply obviously for the other two directions, so in general we can write 



dxj dxi ' 



where Ji. 



dxi 
dxj 



(75) 



(76) 



is the jacobian of this coordinate transformation. This method allows one to have high resolution for 
example near a central object, without however having high resolution anywhere else far away from the 
central object. This is useful in connection with outflows from jets. 

We discuss here one particular application that is relevant for simulating flows in a sphere. It is 
possible to transform a cartesian mesh to cover a sphere without a coordinate singularity. It will turn 
out, however, that there is a discontinuity in the jacobian. We discuss this here in 2-D. We denote the 
coordinate mesh by a tilde, so (i, y) are the coordinates in a uniform cartesian mesh. We want to stretch 
the mesh such that points on the x and y axes are not affected, and that the distance of points on the 
diagonal is reduced by a factor 1/V2 (or by l/\/3 in 3-D). This can be accomplished by introducing new 
coordinates {x, y) as 

'x\ _ fS:^ (i" +y")i/'^ 



where n is a large even number. In the limit n 



yj (i2+y2)l/2' 

oo we may substitute 



(77) 



(i" + y")^/" 



max(|i|, |y|). 

Examples of the resulting meshes for two different values of n are given in figure 
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Figure 11: Examples of the resulting meshes for rt = 8 and n — > oo in which case eq. ( [78|) is used. 

In order to obtain the jacobian of this transformation, dxi/dxj, we have to consider separately the 
cases X > y and x < y. The derivation is given in Appendix and is most concisely expressed in terms 
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of the logarithmic derivative, so 
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if |5| < \y\ 



(79) 



(80) 



where = x 



y^. Note that the jacobian is discontinuous on the diagonals. This is a somewhat 



unfortunate feature of this transformation. It is not too surprising however that something like this 
happens, because the diagonals are the locations where a rotating flow must turn direction by 90° in the 
coordinate mesh. Nevertheless, it is possible to obtain reasonably well behaved solutions; see figure ^ 
for an advection experiment using a prescribed differentially rotating flow. 



— 1 1 r- 



—I 1 1 — 



'./I'll ' I 





Figure 12: Example of an advection experiment on a n = 8 mesh. 

The fluid equations are still solved in rectangular cartesian coordinates, so for example the equation 
Ds/Dt = is solved in the form 

ds ds ds 

^ = -".^--.3^, (81) 

where the spatial derivatives are evaluated according to eq. (75). For the velocity fleld, stress- free bound- 
ary conditions, for example, would be written in the form 

fjUj = 0, inSijfj = 0, (82) 

where Sij is the rate of strain tensor, fj and (pi are the cartesian components {i — x, y, z) of the radial 
and azimuthal unit vectors, i.e. 
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are unit vectors in the r and (f> directions and r — -^/x^ + is the distance from the rotation axis. The 
stress-free boundary conditions are then 

+ yUy — (84) 

and 

(x^ - y^){ux.y + Uy^x) - 2xy{ux,x - Uy^y) = 0. (85) 



4.2 Lagrangian meshes 

We now consider a simple one-dimensional lagrangian mesh problem. Assume that I labels the particle, 
then the lagrangian derivative is 



Ds /9s\ f ds\ f dx\ ds 



Dt \dt) V dt J V dt J dx 



Now, because 

' dx\ Dx 



(86) 



(87) 



we have the well-known equation 

Ds ds ds 

Dt " m'^'^d^' 

As an example we now consider the Burgers equation, 

Du ^d'^u 
D-t=^dx^- 

We now take u{x, t) — u{l{x), t) to be a function of the coordinate variable £ which, in turn, is a function 
of X. The a;-derivatives are obtained using the chain rule, i.e. 

du d£ du u' 

dx dx di x' ' 

and likewise for the second derivative 

d'^u u"x' — u'x" ^^^^ 
dx^ x''^ 



Thus, the Burgers equation can then be written as 

du u"x' — u'x" 

where the x variable is given by 

dx , , 

-=u. (93) 



A solution of these two equations is given in figure 

In the test problem above the initial meshpoint distribution was uniform. Although this is not quite 
suitable for this problem, it shows that subsequently the mesh spacing became narrower still, which 
means that the timestep in now governed by viscosity, 6t < 0.06(5x^;,j/i^, where the numerical factor is 
empirical. However, the mesh spacing does not need to be governed by eq. (p3|), so it is quite possible to 
come up with other prescriptions for the mesh spacing. 
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Figure 13: Solution of the Burgers equation using a lagrangian mesh combined with a sixth order exphcit 
scheme. The values of the mesh Reynolds number vary between 5xU/v = 0.5 and 2.0, where 5x refers 
to the initially uniform mesh spacing. 

Consider as another example the isothermal eulerian equations 

Dm ^ ( d'^u d\npdu\ 



D In p du 
Dt ~ ~dx' 



(95) 



In lagrangian form they take the form 



dt ^ a;' V x'^ x' x' ' ^ 



^ = -- (97) 
%=u. (98) 



The example above demonstrates clearly the problem that lagrangian mesh points can continue to 

pile up near convergence points of the flow. This is a general problem with fully lagrangian schemes. One 
possible alternative is to use lagrangian- eulerian schemes (e.g., Benson 1992, Peterkin, Frese, & Sovinec 
1998, Arber et al. 2001), which combine the advantages of lagrangian and eulerian codes, but involve 
obviously some kind of interpolation. Another alternative is to use a semi-lagrangian code which advects 
the mesh points not with the actual gas velocity it, but with a more independent mesh velocity U. 
Clearly, we want to avoid too small distances between neighboring points, so one could artificially lower 
the effective mesh velocity by involving for example the modulus of the jacobian, |J|, which becomes large 
when the concentration of mesh points is high. Thus, one could choose for example t/ = u/(l + |J|). In 
the present case, |J| = |a;'|~^. In the following we discuss the formalism that needs to be invoked in order 
to calculate first and second derivatives on an advected mesh. 
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4.3 Non-lagrangian mesh advection 

The main advantage of a lagrangian mesh is that it aUows higher resolution locally. Another advantage, 
which is however less crucial, is that the nonlinear advection term drops out. The main disadvantage is 
however that a lagrangian mesh may become too distorted and overconcentrated, as seen in the previous 
section. In this subsection we address the possibility of advecting the mesh with a velocity U that can 
be different from the fluid velocity. This way one can remove the swirl of the mesh by taking a velocity 
that is the gradient of some other quantity, i.e. 



U 



mesh 5 



(99) 



where $mcsh should be large in those regions where many points are needed. One possible criterion would 
be to require that the number of scale heights per meshpoint, \Sx ■ Vlnp|, does not exceed an empirical 
value of 1/3, say. Thus, 3\Sx • Vlnpj < 1 would be a necessary condition. Another possibility would be 
to let $mcsh evolve itself according to some suitable advection-diffusion equation. However, no generally 
satisfactory method seems to be available as yet. In order to calculate the jacobian for the coordinate 
transformation one can make use of the fact that the mesh evolves only gradually from one timestep to 
the next. For a more extended discussion of mesh advection schemes we refer to the article by Dorfi in 
the book by LeVeque et al. (1998). 



4.3.1 Calculating the jacobian 

Initially, at t = 0, we have x ^ x. After the nth timestep, at i = n5t, we calculate the new ai-mesh, 
from the previous one, a;*^"', i.e. 



Sn+i) ^ a;("' +J7(a;("),t) 5t. 



(100) 



Here, a;'-*'' = a; is just the original coordinate mesh. Differentiating the ith component with respect to 
the jth component, as we have done in § 4.1, we obtain 



dx 



(«+i) 



dx 



in) 



dU, 



(n) 



dx 



(n+l) 



dx'j^^^^ dxj 



(n+l) 



dt. 



(101) 



differentiate with respect to the new mesh a;'"+^'. This can be fixed by another factor dx^^^^ /d. 



Thus, we have 



dxP dxi"'> 9C/.(a;("),t) 



dx 



(n+l) 



dx 



(n+l) 



J 3 

This can be written in matrix form, 
where 



dx 



in) 



St ^ 



(n) 



Srk + St 



dx 



(n) 



dx 



in) 



dx 



(n+l) ■ 



where U^"' = [/, (a;("', i). In the expression above we have Ui on the mesh a;("\ but we need to 

,("+!) 

"J 

(102) 

(103) 
(104) 

(105) 



dU,(x'^"'\t) ^ 



dx 



(n) 



is a transformation matrix and 



,(«) 



dx 



(n) 



dx 



(n+l) 



is the incremental jacobian, so J*^"^ = M ^. To obtain the jacobian at t = 2St, for example, we calculate 



dx 



(0) 



dx 



(2) 



dx'',^^ dx^P 



(106) 
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The jacobian at t — nSt is then obtained by successive matrix multiphcation from the right, so 



(O^n+l) ^ ,(0-*n) .(n) 
-'ij •'ik •'kj ' 



(107) 



where and are the full (as opposed to incremental) jacobians at the new and previous 

timesteps, respectively. 

4.3.2 Calculating the second order jacobian 

A corresponding calculation (see Appendix P) for the second derivatives of a function / shows that 



dxidxj dxpdXf 



df 



K 



kij, 



(108) 



where 



d^xk 



is the second order jacobian. Like for the first derivative the second order jacobian can be obtained by 
successive tensor multiplication, 



K 



(O^n+l) 
kij 



uiO-^n) |(n) |(n) |(0^n)^(n) 
'^kpq -"pi -'qj ' -'kl '^lij ' 



(110) 



where K 
and 



(O-^n+l) 
kij 



and k[,° are the second order jacobians at the new and previous timesteps, respectively, 



«■(") 
'^kij 



d'^'-k^ 



dx. 



(Ill) 



is the incremental second order jacobian, which is calculated at each timestep as 



where M was defined in eq. (104) and 



Ui 



l,pq 



d^Ui (a;("),<) 
dxp^^dxlp^ 



(112) 



(113) 



is the second order velocity gradient matrix on the physical mesh. Since M ^ = J*-"' is just the incremental 
jacobian, we can write (112) as 

(114) 



^ki] 



|(") r-i TT ]{■'•■) \{ 
'•'kl ^l,pq-'pi •'qj 



.(n) An) 



Since the expressions ( 110 ) and (114) involve both a multiplication with Jp"''j^"'', we can simplify eq. ( |llO| 



to give 



K 



(O^n+l) 
kij 



i<-(0^n) _ |(0^n) 7-r 
'^kpq •'kl ■^Im"'' ^™.,pq 



An) An) 

pi qj ■ 



(115) 



Here the expression Jfe*;^"^ J;™ is of course the new jacobian, J^^^'^^^^ 
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So in summary, the new first and second order jacobians are obtained from the previous ones via 
the formulae 



|(0^n+l) 



.(O^n) .(n) 
•'ik -'kj ' 



.(O^n+1) 
^kij 



|(") |(") 

pi qj ' 



dxidxi 



dx„dx„ "' "''^ dxk 
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(118) 
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where J = j(o^n+i) ^^^^ ^ |^(o^n+i) j^^g been assumed. 



Since now the mesh is moving in time with the local speed U which is different from the gas velocity 
u, the lagrangian derivative is 

^ = | + («-t/)-V. (120) 

In all other respects the basic equations, written in cartesian form, are still unchanged, provided all 
X, y, and z derivatives (first and second) are evaluated, as in ( |75| ) and ( |108| ), with the components of 
the jacobian. As an example we show in figure ^ the result of a kinematic collapse calculation where 
Dtt/Di = — and Ds/Di = with a smoothed but localized gravitational potential (j>. In figure |l^ we 
compare the results of an eulerian and a lagrangian calculation using the same number of meshpoints. 
Already after some short time the eulerian calculation begins to become underresolved and develops 
wiggles while the lagrangian calculation proceeds without problems. 



4.4 Unstructured meshes 

We now discuss how we can calculate spatial derivatives of our variables from a nonuniformly spaced 
ensemble of points. Consider the function f{x,y,z), which stand for one of the components of a vector 
(velocity or magnetic vector potential) or a scalar, such as In p. We approximate the function f{x,y,z) 
in the neighborhood of the point Xi ~ {xi, yi, Zi) by a multidimensional polynomial of degree N, 

f{x,y,z)^f{xi,yi,z,)+ ^ C4„in{x - x^Yiy - yi)"'{z - z^y, (121) 

l+m+n<N 

where I, m, and n are non-negative integers and c/„m are coefficients that are to be determined separately 
for each point by applying eq. ( |121| ) to all neighboring points Xj. Note that cooo = and does not need 
to be considered. Thus, for each point j we have a system of equations 

^ X! ^i^'nl'^'""^*^ : (122) 

l+m+n<N 

where fij — f{xi,yi,Zi) — f{xj,yj,Zj) and Xij ~ Xi — Xj. This system of equations can be written in 
matrix form 

/a - (123) 

where 1 < (a, /9) < M and M is the spatial dimension of the matrix, which is related to N and the 
dimension as follows: 

( N in 1-D, 

M=< {N +1){N + 2)/2-l in 2-D, (124) 
[ (7V + l)(3iV/2) in 3-D. 
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Figure 14: Example of a collapse calculation. The second row shows only the inner parts with |a;|, \y\ < 
0.01 at the same times as in the upper row. 



When N = 2 the matrix M is given by 
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and 



C = (cioo, CoiO, Cool, C200j Clio, C0205 Coil, C002, Cloi) 



(125) 



(126) 



Here, j„ {n = 1,2,..., M) are the M nearest neighbors of the point i. In general the matrix can be written 
in the form 



M 



(J) 



--a/yaj"'^aj'^ (127) 

where J is the index of the point at which the derivative is to be calculated. The set of exponents Z(/3), 
m(/3), and n(/3) is given here for the case A'' = 4 in 3-D: 

/(/3) = (1,0,0|2,1,0,0,0,1|3,2,1, 0,0,0,0,1, 2|4,3,2,1, 0,0,0,0,0,1, 2,3), (128) 

m{(3) = (0,1, 0|0,1, 2,1, 0,0|0,1, 2,3,2,1, 0,0,0|0,1, 2, 3, 4,3,2, 1,0,0,0,0), (129) 

n{0) = (0,0, 1|0, 0,0, 1,2, 1|0, 0,0, 0,1, 2,3,2, 1|0, 0,0, 0,0, 1,2, 3, 4, 3, 2,1), (130) 

where the vertical bars separate the sets of exponents that correspond to increasing orders. Once the C 
vector has been obtained, the first derivatives of / are simply given by 

cooi- 



df df 

= cioo, 
ox oy 



COlO, 



(131) 



25 




Figure 15: Comparison of lagrangian (+ signs) and eulerian (dots) calculations in the first two plots, and 
later development (last two plots) where the eulerian no longer works. Note that already in the second 
plot the eulerian calculation has developed noticeable wiggles which the lagrangian proceeds without 
problems. 



Likewise, the second derivatives are given by 



d^f d^f d^f _ ,-,oo^ 

C2OO, -TT^ — C02O, TT^ — C002, (132) 



OX'' oy^ az^ 

and the mixed second derivatives are given by 



d^f d^f d^f ,,oo^ 

Clio, TT^^^oii, ■7r^=cioi- (133) 



axay oyoz ozox 

Although this method can be used for meshes that are static in time, it can also be used in connection 
with multi-dimensional lagrangian schemes. In that case there may arise the problem that neighboring 
points get very close together, and so small errors strongly affect the coefficients. A good way out of this 
is to use a few more points and to solve the linear matrix equation using singular value decomposition. 
An example of such a calculation is shown in figure [10, where a passive scalar, with the initial distribution 
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A{x,0) = X, is advected by the velocity, u = r, which in turn is obtained by solving Kepler's equation, 
r = —GMr/r^, using the normalization GM = 1. This windup problem corresponds to the windup of 
initially horizontal magnetic field lines. 



/ - 0.1 ( - 0.5 
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Figure 16: Two-dimensional advection problem on an unstructured lagrangian mesh. The dots indicate 
the 1000 lagrangian particles which constitute the unstructured mesh. 

In diffusivity used in figure |l6| was 77 = 0.02, but due to the coarse resolution and the implicit 
smoothing resulting from the singular value decomposition technique the effective diffusivity is somewhat 
larger. 



5 Implementing magnetic fields 

As mentioned in §^ implementing magnetic fields is relatively straightforward. On the one hand, the 
magnetic field causes a Lorentz force, J x B, where B is the flux density, J = V x B//io is the current 
density, and no is the vacuum permeability. Note, however, that J x B is the force per unit volume, so 
in eq. (||) we need to add the term J x B/p on the right hand side. On the other hand, B itself evolves 
according to the Faraday equation, 

— = -VxE (134) 

where the electric field E can be expressed in terms of J using Ohm's law in the laboratory frame, 
E = —u X B -\- J /a, where a = {r]^o)~^ is the electric conductivity and rj is the magnetic diffusivity. 

In addition we have to satisfy the condition V • B = 0. This is most easily done by solving not for 
B, but instead for the magnetic vector potential A, where B = 'V x A. The evolution of A is governed 



by the uncurled form of eq. (134) 



^^-E-Vc^ (135) 

where is the electrostatic potential, which takes the role of an integration constant which does not 
affect the evolution of B. The choice = is most advantageous on numerical grounds. (By contrast, 
the Coulomb gauge V • j4 = 0, which is very popular in analytic considerations, would obviously be of no 
advantages, since one still has the problem of solving a the solenoidality condition. 

Solving for A instead of B has significant advantages, even though this involves taking another 
derivative. However, the total number of derivatives taken in the code is essentially the same. In fact. 
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when centered finite differences are employed, Alfven waves are better resolved when A is used, because 
then the system of equations for one-dimensional Alfven waves in the presence of a uniform B^q field in 
a medium of constant density po reduces to 

Uz = -^B^oA'' Ay = B^ouz, (136) 

where a second derivative is taken only once (primes denote a;-derivatives). If, instead, one solves for the 
Bz field, one has 

Uz - -^BmB'„ Bz - 5,0^;, (137) 

where a first derivative is applied twice, which is far less accurate at small scales if a centered finite 
difference scheme is used. At the Nyquist frequency, for example, the first derivative is zero and applying 
an additional first derivative gives still zero. By contrast, taking a second derivative once gives of course 
not zero. The use of a staggered mesh circumvents this difficulty. However, such an approach introduces 
additional complications which hamper the ease with which the code can be adapted to other problems. 

Another advantage of using A is that it is straightforward to evaluate the magnetic helicity, {A - B), 
which is a particularly important quantity to monitor in connection with dynamo and reconnection 
problems. 

The main advantage of solving for A is of course that one does not need to worry about the solenoidal- 
ity of the S-field, even though one may want to employ irregular meshes or complicated boundary 
conditions. 

As we have emphasized before, when centered meshes are used, it is usually a good idea to avoid 
taking first derivatives of the same variable twice, because it is more accurate to take instead a second 
derivative only once. For this reason we calculate the current not as J = /Iq x (V x A), but as 

J = li-^[-V^A + V{V ■ A)]. (138) 

Taking the gradient of V • v4 involves of course also taking first derivatives of the same variable twice, 
but these contributions are canceled by corresponding components of the V^v4 term. There are some 
advantages relying here on the numerical cancellation, which is of course not exact. The reason is that the 
full V^A term is important when used in the magnetic diffusion term. If the diagonal terms, / dx^ , 
d^Ay/dy^, and d^Az/dz^, which would all drop out analytically, were taken out there would be no 
diffusion of A in the direction of A. 

There is one more aspect that is often useful keeping in mind. There is a particular gauge that allows 
one to rewrite the uncurled induction equation in such a form that the evolution of A is controlled by 
the advective derivative of A. The calculation is easy. Write the induction term u x B va. component 
form and express B in terms oi A, so 

(m X V X A)i — Uj{diAj — d-jAi) — di{ujA.j) — Ajd^Uj — UjdjAi. (139) 

Here the last term contributes to the advective derivative, the first term can be removed by a gauge 
transformation and the middle term is a modified stretching term, so the induction equation takes the 
form 

DA 

— i = -AjdiUj - 77/io Jj. (140) 
This gauge was used by Brandenburg et al. (1995) in order to treat a linear velocity shear using pseudo- 



periodic (shearing box) boundary conditions. The formulation (14C) can also be useful when solving the 
induction equation using lagrangian methods. Note, however, that the nonresistive evolution of A differs 
from that of B in that the indices of the matrix Uy = dui/dxj are interchanged and that the sign is 
different; positive for the -B-equation, 

Di?- 

—1 +UijBj + other terms, (141) 
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and negative for the A-equation, 

DA; 



= -Ajiiji + other terms. (142) 

These two formulations are particularly advantageous when the velocity has a constant gradient, as in 
the case of linear shear. In local simulations of accretion discs, for example, the shear component is 

DA^ 3, 



Uy{x) = —^flx, SO \Jyx — — fi^, and all other Uy vanish. Hence 



flAy + other terms (143) 



Dt 2 
for the j4-formulation, or 

T)]3 3 

y- = — flB^ + other terms (144) 

Dt 2 ' 

for the B-formulation. In these two formulations all dependent variables are clearly periodic (or rather 

pseudo-periodic), so there is no term that is explicitly non-periodic such as Uy{x) = — |f2a;. In the 

following, whenever magnetic fields are present, we use the A-formulation, mainly because it guarantees 

the solenoidality of B everywhere (including the boundaries), and also because it is easy to use. 



Cache-efRcient coding 

Unlike the CRAY computers that dominated supercomputing in the SOties and early 90ties, all modern 
computers have a cache that constitutes a significant bottleneck for many codes. This is the case if large 
three-dimensional arrays are constantly used within each time step. The advantage of this way of coding 
is clearly the conceptual simplicity of the code. A more cache-efficient way of coding is to calculate 
an entire timestep (or a corresponding substep in a three-stage 2N Runge-Kutta scheme) only along a 
one-dimensional pencil of data within the box. On Linux and Irix architectures, for example, this leads 
to a speed-up by 60%. An additional advantage is a drastic reduction in temporary storage that is needed 
for auxiliary variables within each time step. 



6 Application to astrophysical outflows 
6.1 The isothermal Parker wind 

Before discussing outflows from accretion discs it is illuminating to consider first the one-dimensional 
example of pressure-driven outflows in spherical geometry. A particularly simple case is the isothermal 
wind problem, which is governed by the equations 

at or or Or 

where Cs is the isothermal sound speed (assumed constant), M is the mass loss rate, and ^(r) is a 
prescribed function of position, normalized such that J47rr^^(r)dr — 1, and non- vanishing only near 
r = 0. For a point mass the gravity potential $ would be —GM/r, but this becomes singular at the 
origin. Therefore we use the expression $ = -GM/{r'^ + rjf)i/" instead, where we choose n = 5 in all 
cases, and 1/ro gives the depth of the potential well. In figure |l^ we show radial velocity and density 
profiles for different values of M . Note that the velocity profile is independent of the value of Af , but the 
density profile changes by a constant factor. In the steady case the equations can be combined to 

dl^^2c^_GM^ (147) 

so the sonic point, |u| — Cg, is at r = = GM/2Cg. In figure |lj we have chosen GM — 2 and Cg — 1, so 
= 1, which is consistent with the graph of u. 
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Figure 17: Isothermal Parker wind solutions for different values of M . 
independent of the value of M . GM = 2, Cg = 1, ro = 0.4. 



Note that the u profile is 



6.2 The poly tropic or adiabatic wind 

In the following we make the assumption that the entropy is constant. In that case it is particularly 
useful to solve for the potential enthalpy^ H = h + which varies much less than either h or ^. Using 
H as dependent variable is particularly useful if one solves the equations all the way to the origin, r — 0, 
where $ tends to become singular (or at least strongly negative if a smoothed potential is used) . In terms 
of H the governing equations are 

du du dH ^ 

uM^, (148) 



dt 



dr 



dH 
'dt 



dH 
or 



or 



dr 
1 d 



(r^u) 



M^(r) 



(149) 



where H — h + ^ is the potential enthalpy, h — p/p + e is the enthalpy, and Cg ~ (7— l)h — (7— l)(i/ — $) 
for a perfect gas, where Cg is the adiabatic sound speed and h ~ CpT is the enthalpy. These equations 
are also valid in the nonisothermal case (7 7^ 1). The isothermal case may be recovered by putting 7=1 
and replacing h by Cg In p. In figure |l^ we show solutions for different values of M and 7 = 5/3. Again 
we put GM = 2 and Cgo = 1. 

We note that, depending on the strength of the mass source, the polytropic wind problem allows a 
variety of different velocity and Mach number profiles, whereas for the isothermal wind problem there 
was only one solution possible, independent of the strength of the mass source. The velocity profile was 
always the same and also the density was the same up to some scaling factor that changes with M. This is 
connected with the additional degree of freedom introduced through the polytropic constant K = p/rho"' . 
Since Cg is no longer constant, the position of the sonic point is no longer fixed and different solutions are 
possible. 



In figure 19 we show solutions where M is kept constant, but the depth of the potential well, GM/rQ, 
is changed by varying the value of tq. Note that the deeper the potential well, the higher the wind speed. 
The density far away from the source is then correspondingly smaller, so as to maintain the same mass 
flux. 

As we have seen in § 3.6, a polytropic equation of state is unphysical. Therefore we now consider the 
case where the energy equation is included. To be somewhat more general we consider first the basic 
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Figure 18: Polytropic Parker wind solutions for different values of M . GM — 2, Cso = 1, ''o = 0.4. 
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Figure 19: Density p, velocity u, and Mach number Ma = u/cg for the polytropic Parker wind solutions 
for different values of tq. GM = 2, Cso = 1, M = 10. 



equations in conservative form with mass, momentum and energy sources included, i.e. 



dt 



dx-i 



dt 



\pu^ + pe) + ^— [uj {\pu'^ + ph) - UiTij] = i?^. 



(150) 

(151) 
(152) 



where Af, J, and E are the rates of mass, momentum and energy injection into the system, — 2vpSij 
is the viscous stress tensor, and Sij is the (traceless) rate of strain tensor; see eq. (p^). Rewriting the 
energy equation in nonconservative form we have 



2£ + Ev.«^Tgi = 2.S2 

Dt p m 



E-u-[l- uM^ - (iw^ + e) 



M 



Or). 
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P 



(153) 
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which can also be rewritten in terms of entropy, so the final system of nonconservative equations with 
source terms is 



Dlnp 
Bt 



V ■u = M 



I - uM 



E-ui+ {\v? - e) M 



p 



(154) 
(155) 
(156) 



where T can be replaced by Cg/{'j — 1) (remember that Cp = 1), and Cg is given by (0). 




Figure 20: Wind solutions for different values of E and M — 10. Note that the solution with E ~ 35 
is quantitatively very similar to the polytropic solution with the same value of M. GM = 2, Cso — 1, 
ro = 0.4. 



In figure 20 we present solutions of Eqs. (154)-( |156D for different values of E and M. The main effect 
of varying the value of E is to change the value of the entropy in the wind. Outside the acceleration 
region, however, the value of the entropy is fairly constant, so the polytropic assumption appears to be 
reasonably good here. 

While outflows of some very early-type stars are driven mostly by the I term (resulting from the 
radiation pressure in lines), the winds of cool stars are driven mostly by the E term (resulting from the 
hot coronae). Similar differences may also explain why some jets are massive (stellar jets, for example), 
whilst others are not (jets from active galactic nuclei, for example, or those anticipated in gamma-ray 
bursters). 

6.3 Relevance to outflows and jets 

The pressure-driven outflows discussed in the previous section may take the form of more collimated 
outflows once a magnetic field is involved. This applies to the case of magnetized accretion discs. These 
discs are generally magnetized both because of dynamo action within the disc and because of external 
fields that were dragged into the disc from outside due to the accretion flow. 

At least in some types of jets the outflows may be driven by hot coronae. Other possibilities for 
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driving outflows involve the magneto-centrifugal effect. It is well-known that outflows can be driven 
from a magnetized disc if the angle between the field and the disc is less than 60° (Blandford & Payne 
1982). Recent work in this field was directed to the question whether this angle is the result of some 
self-regulating process (Ouyed, Pudritz, & Stone 1997, Ouyed & Pudritz 1997a, b, 1999) and whether it 
can be obtained automatically from a dynamo operating within the disc (Campbell 1999, 2000, Dobler 
et al. 1999, M. v. Rekowski, Riidiger, & Elstner 2000). This latter question is particularly interesting 
in view of the fact that jets in star-forming regions are not really pointing in a similar direction (e.g. 
Hodapp & Ladd 1995), as one might expect from jet models that start off with a prescribed large scale 
field. 




Figure 21: Poloidal velocity vectors and poloidal magnetic field lines superimposed on a grey-scale image 
of the logarithmic density. Dark means high density. The thick dashed line denotes the location where 
the poloidal flow speed equals the fast magnetosonic speed. The thin solid line gives the location of the 
disc surface. The slight asymmetry in the field is a relic from the mixed-parity initial condition. [Adapted 
from Brandenburg (2000).] 

In figure |l] we present a particular model of Doblcr et al. (1999) and Brandenburg (2000); see Bran- 
denburg et al. (2000) for a full account of this work. In these models the outfiow is driven by mass 
sources whose strength is proportional to the local density deficit relative to that in the original equilib- 
rium solution of the disc. Such a density deficit was initially caused by slow gas motions that resulted 
from an instability of the initial equilibrium solution, because a cool disc embedded in a hot corona is 
nonrotating outside the disc, and it is the resulting vertical shear profile that causes the instability (cf. 
Urpin & Brandenburg 1998). At later times, of course, the outflow makes the corona corotating, but by 
that time the outflow is driven by a persistent density deficit in the disc relative to the initial references 
solution. 

In this model the magnetic field was generated by an a— $7 dynamo operating within the disc. However 
a is negative in the upper disc plane (see Brandenburg et al. 1995), and then the most preferred field 
geometry is dipolar (Campbell 1999, v. Rekowski, Riidiger, & Elstner 2000). The field parity is sensitive 
to details in the disc physics assumed in the particular model (aspect ratio, disc thickness, the presence of 
outflows, and the conductivity in the disc and the exterior). Nevertheless, both dipolar and quadrupolar 
fields are equally well able to contribute to wind launching, at least in the outer parts of the disc where the 
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angle between the field and the disc plane is less than 60° , the critical angle for magneto-centrifugal wind 
launching (Blandford & Payne 1982). We note, however, that the more detailed analysis of Campbell 
(1999) suggests that the critical angle can be significantly smaller. 

In our models the outflow is only weakly coUimated (if at all). This is probably connected with the 
fact that here the fast magnetosonic surface is rather close to the disc surface, making it difhcult for the 
field to become strong enough to channel the magnetic field. Instead, the field lines themselves are still 
being controlled too strongly by the outflow. However, outflows with rather large opening angles are 
actually seen in some star-forming regions; see Greenhill et al. (1998). 

While most of the disc mass is ejected in a cone of half-opening angle around 25°, most of the disc 
angular momentum is ejected at rather low latitudes, almost in the direction of the disc plane away from 
the central object. The timescales for these various processes are comparable. In flgure ^we show the 
azimuthally integrated mass flux, angular momentum flux, and magnetic (Poynting) flux as a function of 
polar angle, and compare with a nonmagnetic run. We flnd that in the magnetic run the outflow is more 
strongly concentrated towards the axis. Also, the amount of angular momentum loss (dash-dotted line) 
is larger when the disc is magnetized. We emphasize in particular that in the magnetic run significant 
amounts of magnetic field are eject from the system. In the following section we discuss the significance of 
such magnetic flux ejection for magnetizing the interstellar medium into which the outflow is streaming. 
This discussion is similar to a corresponding discussion for the contamination of the intergalactic medium 
via outflows from active galactic nuclei (Brandenburg 2000). 
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Figure 22: Comparison of the angular dependence of azimuthally integrated fluxes for magnetic and 
nonmagnetic outflows. The solid line refers to mass flux, the dashed line to angular momentum flux, and 
the dash-dotted line (in the second panel) corresponds to the Poynting flux. The units of all quantities 
are thus l/[t]. 



6.4 Magnetic contamination from outflows 

It may at first appear somewhat unrealistic to expect significant magnetization of the interstellar medium 
from outfiows. However, the following calculation shows that the effect may be quite significant. Assume 
that every star did undergo a phase of strong accretion with associated outfiows, so iV = 10^^ for the whole 
galaxy. The duration of intense outfiow activity is 10^ years, say, but it could even be 10^ years. The 
magnetic luminosity is Lmag = COSMwCg (Brandenburg et al. 2000), where Cg ~ lOkm/s is the average 
sound speed of the interstellar medium, and = O.lMd (see Pelletier & Pudritz 1992), where ~ 
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10 "^Mq/ yr is a conservative estimate for the disc accretion rate. Again this value may be larger. With 
the above numbers the magnetic luminosity from all N = 10 sources is then iVLmag = 7 x 10'^^ erg/ s 
and the total energy output delivered from all stars at some early point in the life time is therefore 
Emag ~ TNLniag — 2 X 10^^ erg. Diluting this over a volume of a galaxy of 300 kpc'^ (radius 10 kpc, 
height 1 kpc) gives 2 x 10^^^ erg/ cm'^. Multiplying this by Stt and taking the square root gives 0.2/xG. 
Expressed more concisely in a formula we have for the rms magnetic field strength 



where the efficiency factor Fpoy/Jkin (=0.05 in our model) may be lower in systems where the disc 
dynamo is less strong. 

The parameters for a corresponding estimate for outflows from young galactic discs (active galactic 
nuclei) are as follows. Assuming N ~ 10'' galaxies per cluster, each with O.lM0/yr = lO^^g/s, 

and Cs = 1000 km/s for the sound speed in the intracluster gas, the rate of magnetic energy injection for 
all galaxies together is Lmag = 10'*^crg/s. Distributing this over the volume of the cluster of F ^ 1 Mpc'', 
and integrating over a duration of St = 1 Gyr, this corresponds to a mean magnetic energy density 
of (B^/Stt) « 10-13 erg/cm^, so (B2)1/2 ~ iQ-s G, which is indeed of the order of the field strength 
observed in galaxy clusters. We note that our estimate has been rather optimistic in places {M„ could 
be lower, or the relevant 5t could be shorter, for example), but it does show that outflows are bound to 
produce significant magnetization of the intracluster gas and the interstellar medium (see also Volk & 
Atoyan 1999). In the latter case it will provide a good seed field for the galactic dynamo. A dynamo 
is still necessary to shape the magnetic field and to prevent if from decaying in the galactic turbulence. 
Similarly, many galaxy clusters undergo merging and this too can enhance and reorganize the magnetic 
field. The necessity for a recent merger event would also be consistent with the fact that not all halos are 
observed to have strong magnetic fields. Recent simulations by Roettiger, Stone & Burns (1999) suggest 
that after a merger the field strength may increase by a factor of at least 20 (and this value increases 
with improving observational resolution). 

As an alternative consideration for causing the magnetization in clusters of galaxies, primordial mag- 
netic fields are sometimes discussed. There are numerous mechanisms that could generate relatively 
strong fields at an early time, for example during inflation (age ^ IQ-^^ s) or during the electroweak 
phase transition (age ~ 10-^° s). Such fields would now still be at a very small scale if one considers only 
the cosmological expansion. However, depending on the degree of magnetic helicity in this primordial 
field, the magnetic energy can be transferred to larger scales that are now on the scale of galaxies. For a 
recent discussion of these results see Brandenburg (2001a). 

7 Hydromagnetic turbulence and dynamos 

As mentioned in the beginning, accurate high order schemes are essential in all applications to turbulent 
flows. Nevertheless, we should mention that one often attempts solutions of the inviscid and nonresistive 
equations using low-order finite differences combined with monotonicity schemes that result in some 
kind of effective diffusion. The piece-wise parabolic method (PPM) of Colella & Woodward (1984) is an 
example. However, unlike the Smagorinsky scheme (see Chan & Sofia 1986, 1989, Steffen, Ludwig, & Kriifi 
1989, Fox, Theobald. & Sofia 1991 for applications to convection simulations), PPM and similar methods 
cannot be proven to converge to the original Navier-Stokes equation in the limit of infinite resolution. 
Nevertheless, they are rather popular in astrophysical gas simulations. These schemes are rather robust 
and have also been applied to high resolution simulations of compressible turbulence (Porter, Pouquct, 
& Woodward 1992, 1994). While the results from those simulations are generally quite plausible, the 
power spectrum shows a. subrange at large wave numbers, which is still not fully understood. This 
was sometimes regarded as an artifact of PPM, and should therefore only occur at small scales. However, 
as the resolution was increased further (up to 1024'^), the subrange just became more extended. 




(157) 
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A similar feature was found in cascade models of turbulence when the ordinary diffusion operator 
was replaced by a —V'' "hyperdiffusion" operator (Lohse & Miiller-Groeling 1995). Whatever the outcome 
of this puzzle is, it is clear that with schemes that cannot be proven to converge to the actual Navier- 
Stokes equations in the limit of infinite resolution, there would always remain some uncertainty and 
debate. On the other hand, especially in the incompressible case the use of hyperviscosity does generally 
allow the exploration of larger Reynolds numbers and broader inertial ranges. 

MHD simulations with the highest resolution to date have been performed by Biskamp & Miiller 
(1999), who considered decaying turbulence with and without magnetic helicity. They found that in the 
presence of magnetic helicity the magnetic energy decay is significantly slower. In particular, they found 
the magnetic energy decays like t~^/^, as opposed to t^^ found earlier by Mac Low, Klessen, & Burkert 
(1998) for compressible turbulence. 

Before we start discussing dynamo action in turbulence simulations representative of more astrophys- 
ical settings, such as accretion discs and the solar convection zone, let us first illustrate the mechanism of 
the inverse cascade that is believed to be an important ingredient of large scale magnetic field generation. 



7.1 Isotropic MHD turbulence 

Most developments in the theory of turbulence have been carried out under the assumptions of homo- 
geneity and isotropy. This is certainly true of the work on the inverse cascade (or turbulent cascades 
in general), but it is also true of much of the work on the a-effect which - like the inverse cascade - 
describes the generation of large scale fields. However, unlike the inverse cascade process, the energy 
comes here directly from the velocity field at the scales of the energy-carrying eddies and not from the 
velocity and magnetic field at successively smaller scales, which are usually larger than the scale of the 
energy- carrying eddies. 

It is not easy to see whether any of these effects is actually responsible for the large scale field 
generation in astrophysical bodies or even the simulations. In simulations of accretion disc turbulence 
there is certainly some evidence for the presence of an a-effect, but it is extremely noisy (Brandenburg 
et al. 1995, Brandenburg & Donner 1997, Ziegler & Riidiger 2000). Evidence for the inverse magnetic 
cascade comes mostly from the magnetic energy spectra (Balsara & Pouquet 1999, Brandenburg 2001b), 
which show a marked peak at large scales, but this is convincing only in cases where the flow is driven 
at a wavenumber that is clearly larger than the smallest wavenumber in the box. In practice, e.g. in 
convectively driven turbulence, the flow is driven at all scales including the large scale making it difficult 
to see a marked peak at the smallest wavenumber (see a corresponding discussion in Meneguzzi & Pouquet 
1989). 

From the seminal papers of Frisch et al. (1975) and Pouquet, Frisch, Leorat (1976) it is clear that 
amplification of large scale fields can also be explained by an inverse cascade of magnetic helicity. In 
those papers the authors also showed that the inverse cascade is a consequence of the fact that the 
magnetic helicity, (A • B), is conserved by the nonresistive equations. (A is the magnetic vector potential 
giving the magnetic field as B = V x A.) The inverse magnetic cascade effect too is rather difficult to 
isolate in simulations of astrophysical turbulence. However, under somewhat more idealized conditions, 
for example when magnetic energy is injected at high wave numbers, one clearly sees how the magnetic 
energy increases at large scales; see figure Further details of this model have been published in the 
proceedings of the helicity meeting in Boulder (Brandenburg 1999). 

In the model considered above the fiow was forced magnetically. This may be motivated by the 
recent realization that strong magnetic field generation in accretion discs can be facilitated by magnetic 
instabilities, such as the Balbus-Hawley instability. Other examples of magnetic instabilities include 
the magnetic buoyancy instability, which can lead to an a-effect (e.g. Brandenburg & Schmitt 1998, 
Thelen 2000), and the reversed field pinch which also leads to a dynamo effect (e.g. Ji et al. 1996). Before 



returning to the accretion disc dynamo in § 7^9 we should emphasize that strong large scale field generation 
is also possible with purely hydrodynamic forcing. Simulations in this type were considered recently by 
Brandenburg (2001b). There are many similarities compared with the case of magnetic forcing. The 
evolution of magnetic energy spectra in the presence of hydrodynamic forcing is shown in figure UA. Like 
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Figure 23: Spectral magnetic energy, Eyi{k^ t), as a function of wavcnumber k for different times: dotted 
lines are for early times {t — 2,4,10,20), the solid and dashed lines are for t — 40 and 60, respectively, 
and the dotted-dashed lines are for later times {t — 80, 100, 200, 400). Here magnetic energy is injected at 
wavenumber 10. Note the occurrence of a sharp secondary peak of spectral magnetic energy at fc = 10. 
By the time the energy at /c = 1 has reached equipartition the energies in fc = 2 and fc = 4 become 
suppressed. 

in the case of magnetic forcing (figure ^) there are marked peaks both at the forcing scale and at the 
largest scale of the box. Furthermore, the evolution of spectral energy at the largest scales shows similar 
behavior: the magnetic energy with wavenumber k — 8 increases, reaches a maximum, and begins to 
decrease when the magnetic energy at fc = 4 reaches a maximum. The same happens for the next larger 
scales (wavenumbers fc = 4 and 2, until the scale of the box (with fc = 1) is reached. 

The suppression of magnetic energy at intermediate scales, 2 < fc < 8, is quite essential for the devel- 
opment of a well-defined large scale field. In a recent letter Brandenburg & Subramanian (2000) showed 
that this type of self-cleaning effect can also be simulated by using ambipolar diffusion as nonlinearity and 
ignoring the Lorentz force altogether. Without any nonlinearity, however, there would be no interaction 
between different scales and the magnetic energy would increase at all scales, especially at small scales, 
which would soon swamp the large scale field structure with small scale fields. 

The model presented in figure ^ has large scale separation in the sense that there is a large gap 
between the forcing wavenumber (fc = fcf = 30) and the wavenumber of the box (fc = fci = 1). One sees 
that during the growth phase there is a clear secondary maximum at fc = 7. This is indeed expected for 
an dynamo, whose maximum growth rate is at fcmax = ■^<^/vt, where ?7t is the total (turbulent plus 
microscopic) magnetic diffusion coefficient. 

The disadvantage of a high forcing wavenumber is that for modest resolution (here we used 120'^ 
meshpoints) no inertial range can develop. This is different if once forces at fcf = 5, keeping otherwise 
the same resolution. In figure ^ we show spectra for different cases with fcf = 5 where we compare the 
results for different values of the magnetic Reynolds and magnetic Prandtl number. In figure ^ we show 
cross-sections of one field component at different times. In this model (Run 3 of Brandenburg 2001b) the 
forcing is at fcf = 5, so there is now a clear tendency for the build-up of an inertial range in 8 < ^ < 25. 

7.2 The inverse cascade in decaying turbulence 

We now turn to the case of decaying turbulence, which is driven only by an initial kick to the system. 
There are several circumstances in astrophysics where this could be relevant: early universe, neutron 
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Figure 24: Left: magnetic energy spectra for a run with forcing at fc = 30. The times range from 
(dotted line) to 10, 30, 290 (solid lines). The thick solid line gives the final state at t = 1000. Note 
that at early times the spectra peaks at fcmax ~ 7. The and k^'^l'^ slopes are given for orientation as 
dash-dotted lines. Right: evolution of spectral magnetic energy for selected wavenumbers in a simulation 
with hydrodynamical forcing at k = 30. 

stars, and mergers of galaxy clusters. In all those cases one is interested in the development of large 
scale fields. In the context of the early universe the possibility of energy conversion from small to large 
scale fields was pointed out by Brandenburg, Enqvist, & Olesen (1996) who found that fields generated 
at the horizon scale of 3 cm after the electroweak phase transition would now have a scale on the order 
of kiloparsecs, even though the cosmological expansion alone would only lead to scales on the order of 
1 AU. These results were only based on either two-dimensional simulations or three-dimensional cascade 
model calculation (e.g. Biskamp 1994). Therefore we now turn to fully three-dimensional simulations. 

In the absence of any forcing and with no kinetic energy initially an initial magnetic field can only 
decay. However, if initially most of the magnetic energy is in the small scales, there is the possibility 
that magnetic helicity and thereby also magnetic energy is transferred to large scales. This is exactly 
what happens (figure , provided there is initially some net helicity. The inset of figure |2^ shows that 
in the absence of initial net helicity the field at large scales remains unchanged, until diffusion kicks in 
and destroys the remaining field at very late times. 

If the magnetic field has the possibility to tap energy also from the large scale velocity the situation 
is somewhat different again and there is the possibility that a large scale magnetic field can also be 
driven without net helicity. In that case the large scale field can increase due to dynamo action from the 
incoherent a — 17-effect (Vishniac & Brandenburg 1997). In astrophysical settings there is usually large 
scale shear from which energy can be tapped. Before we discuss simulations with imposed shear in more 
detail we first present a simple argument that makes the link between the inverse cascade and helicity 
conservation. 

7.3 The connection with magnetic hehcity conservation 

In the following we give a simple argument due to Frisch et al. (1975) that helps to understand why the 
magnetic helicity conservation property leads to the occurrence of an inverse cascade. We define in the 
following magnetic energy and helicity spectra, M{k) and -ff (fc), respectively. Now, because of Schwartz 
inequality, we have 

|B(fc)|2 ^\{ikxA)-B\> \k\\A ■ B\ (158) 
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Figure 25: Comparison of time averaged magnetic energy spectra for Runs 1-3 {t = 600 — 1000) with 
a non-averaged spectrum for Run 5 (large magnetic Prandtl number) taken ai t — 1600. To compensate 
for different field strengths and to make the spectra overlap at large scales, two of the three spectra 
have been multiplied by a scaling factor. There are clear signs of the gradual development of an inertial 
subrange for wavenumbers larger than the forcing scale. The k~^/^ slope is shown for orientation. The 
dissipative magnetic cutoff wavenumbers, (J^/?7^)^/'', are indicated by arrows at the top. 



we have a lower bound on the spectral magnetic energy at each wavenumber k = k. In terms of shell 
integrated magnetic energy and helicity spectra this corresponds to 

M{k) > ^k\H{k)\, (159) 

where the 1/2-factor comes simply from the 1/2-factor in the definition of the magnetic energy. Assuming 
that two wave numbers q and p interact such that they produce power at a new wave number k, then 

M{p) + M{q) = M{k), H{p) + H{q) = H{k). (160) 

For simplicity we consider the case p = q, so 

2M{p) = M{k), 2H{p) = H{k). (161) 

Assume also that initially the constraint was sharp (maximum helicity), then 

M{p) = \pH{p). (162) 

Now, from the constrain again we have 

\kH{k) < M{k) = 2M{p) = pH{p) = \pH{k), (163) 
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Figure 26: Gray-scale images of cross-sections of Bx{x, y, 0) for Run 3 of Brandenburg 2001b at different 
times showing the gradual build-up of the large scale magnetic field after t = 300. Dark (light) corresponds 
to negative (positive) values. Each image is scaled with respect to its min and max values. 

so 

k < p, (164) 

that is the wave number of the target result must be larger or equal to the wave numbers of the initial 
field. 

The argument given above is of course quite rough, because it ignores for example the detailed 
angular dependence of the wave vectors. This was taken into account properly already in the early paper 
by Pouquet, Prisch, & Leorat (1976), but this approach was based on closure assumptions for the higher 
moments, which is in principle open to criticism. Thus, numerical simulations, like those presented above, 
are necessary for an independent confirmation that the inverse cascade really works. In this connection 
one should mention that there are some parallels with the inverse cascade of enstrophy in two-dimensional 
hydrodynamic (nonmagnetic) turbulence. In that case the enstrophy (i.e. the mean squared vorticity) 
is conserved because of the absence of vortex stretching in two dimensions. The inverse hydrodynamic 
cascade has some significance in meteorology and perhaps in low aspect ratio convection experiments, 
where one finds a peculiar energy and entropy spectrum that is referred to as Bolgiano scaling; see 
Brandenburg (1992) and Suzuki & Toh (1995) for corresponding shell model calculations and Toh & lima 
(2000) for direct simulations. 
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Figure 27: Power spectra of magnetic energy (solid lines) and kinetic energy (dotted lines) in a decay run 
with helicity. The left hand panel is for a case where the flow is only driven by an initial helical magnetic 
field. In the right hand panel the field is weak and governed by strong decaying fluid turbulence. The 
inset shows both velocity and magnetic spectra in the same plot. The Prandtl number v/rj is equal to 
one in both cases, but the mesh Reynolds number, which is kept constant at all times, is different in the 
two cases: 20 in the left hand panel and 50 in the right hand panel. The times are 0, 0.01, 0.1, etc, till 
t — 10^ in the left hand panel and t = 10"^ in the right hand panel. 



7.4 Inverse cascade or a-effect? 



In § 7.1 we made a distinction between inverse cascade and a-effect in the sense that, although both lead 
to large scale field generation, in the inverse cascade there is a gradual transfer of magnetic helicity and 
energy to ever larger scales, whereas the a-effect produces large scale magnetic fields directly from small 
scale fields. Thus, the distinction is really one between local and nonlocal inverse cascades. 



In figure 28 we show the normalized spectral energy transfer function T{k,p,t) for k = 1 and 2 as a 
function of p, and at different times t. The index k signifies the gain or losses of the field at wavenumber 
fc, and the index p indicates the wavenumber of the velocity from which the energy comes from. This 
function shows that most of the energy of the large scale field at fc = 1 comes from velocity and magnetic 
field fluctuations at the forcing scale, which is here fc = fcf = 5. At early times this is also true of the 
energy of the magnetic field at fc = 2, but at late times, t = 1000, the gain from the forcing scale, fc = 5, 
has diminished, and instead there is now a net loss of energy into the next larger scale, fc = 3, suggestive 
of a direct cascade operating at fc = 2, and similarly at fc = 3. 

Based on these results we may conclude that in the saturated state the magnetic energy at fc = 1 
is sustained by a nonlocal inverse cascade from the forcing scale directly to the largest scale of the box. 
This is characteristic of the a-effect of mean-field electrodynamics, except that here nonlinearity plays 
an essential role in isolating the large scale from the small scale 'magnetic trash', as Parker used to say. 

A closer look at figure ^ where fc fcf 30, suggests that once the scale separation is large enough 
the energy is at first transferred not to the scale of the box, but instead to a somewhat smaller scale (here 
at wavenumber fc = 7). Following the corresponding discussion in Brandenburg (2001b), this wavenumber 
is close to the wavenumber, fcmax = ||a|/?7T, where the a^ dynamo grows fastest. 

In the following section we address the issue of magnetic helicity conservation which has important 
consequences for the timescale after which the large scale field begins to develop. This has also a bearing 
on the widely discussed controversy of the so-called 'catastrophic a-quenching' of Vainshtein & Cattaneo 
(1992). 
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Figure 28: Spectral energy transfer function T{k,p,t), normalized by (B^) for three different times, for 
fc = 1 and 2. Run 3 of Brandenburg (2001b). 



7.5 Approximate helicity conservation 

The magnetic helicity, H ^ {A ■ B), is conserved by the nonresistive MHD equations. For a closed or 
periodic box (A • B) satisfies the equation 

^^{A-B) = -2Tjfio{J-B), (165) 

■where (J-B) is the current helicity, and angular brackets denote volume averages. Note that for a periodic 
box {A ■ B) is gauge invariant, i.e. {A ■ B) does not change after a gauge transformation, A — > A + Wtp. 
This is a direct consequence of the solenoidality of the magnetic field, because ( Viy9 ■ B) = — {ip'V ■ B) = 
owing to V • S = 0. 

In order to judge whether {A ■ B) is small or large we calculate the length scale 

in^\{A-B)\/{B^). (166) 

In figure ^ we see that the evolution of proceeds in three distinct phases: (i) a very short period 
(t < 1) where is very small and comparable to the numerical noise level, so magnetic helicity almost 
perfectly conserved, (ii) an intermediate interval (2 < < < 200) where t-a is much larger, but still only 
roughly equal to the mesh size of the calculation, and then (iii) a regime where is of order unity. 
The latter is only possible because of the presence of helicity in the system, which leads to a large scale 
magnetic field configuration that is nearly force-free. 



7.6 Resistively limited growth of the large scale field 

The approximate conservation of magnetic helicity has an important consequence for the generation of 
large scale fields: in order to build up a large scale field with magnetic helicity one has to change the value 
of {A ■ B) from its initial value of zero to its certain final value. This final value of {A ■ B) is such that the 
length scale €h is close to the maximum value possible for a certain geometry. It also implies that when 
interpreting the results in terms of mean- field electrodynamics (a-effect and turbulent diffusion), the a- 
effect must be quenched early on, well before the field has reached its final value. Nevertheless, this strong 
quenching, which was first anticipated by Vainshtein & Cattaneo (1992) and later confirmed numerically 
by Cattaneo & Hughes (1996), does not prevent the field from reaching its final super-equipartition field 
strength. 
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Figure 29: Evolution of the (negative) magnetic helicity length scale in a double-logarithmic plot. 
Note the presence of three distinct phases: very approximate helicity conservation near zero, followed 
by a phase of larger magnetic helicity scale (three orders of magnitude), and finally a phase where the 
magnetic helicity scale has reached the scale of the box. 



Before we come to the details we mention already now that the helicity constraint is probably too 
severe to be acceptable for astrophysical conditions, so one must look for possible escape routes. The most 
plausible way of relaxing the helicity constraint is in allowing for open boundary conditions (Blackman 
& Field 2000, Kleeorin et al. 2000), but the situation can still be regarded as inconclusive. Given that 
much of the work on large scale dynamos so far assumes periodic boundaries, we shall now consider this 
particular case in more detail. In the periodic case the final field geometry can be, for example, of the form 
B = _Bo(cosfciz,sinfciz,0), where fci = 1 is the smallest possible wavenumber in the box, i?o is the field 
amplitude, and {B ) = Bq. Alternatively, the field may vary in the x ot y direction, and there may be an 
arbitrary phase shift; examples of these possibilities have been reported in Brandenburg (2001b). Anyway, 
for B = Bo{coskiz,sinkiz,0) the corresponding vector potential is A = — (Bo/A;i)(cos fciz, sinfciz, 0) + 
V0, where is an arbitrary gauge which does not affect the value of {A ■ B). In this example we have 

{A.B) = -{B^)/h, (167) 



where we have included the fci factor, even though in the present case ki — l. (The minus sign in eq. ( |167| ) 
would turn into a plus if the forcing had negative helicity.) The mean current density is given by 
J = —Boki(coskiz,shikiz,0), so the current helicity of the mean field is given by 

/io(J-B) =-fci(:B'). (168) 



Before we can use eqs. (167) and (168) in eq. (165) we need to relate the magnetic and current helicities of 



the mean field to those of the actual field. We can generally split up the two helicities into contributions 
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from large and small scales, i.e. 



(A- B) = {A- B) + (a- b), 
{J-B)^(j-B) + {j-b). 



(169) 
(170) 



As t he la rge scale magnetic field begins to saturate, the magnetic helicity has to become constant and so 
eq. (165) dictates that (J • B) must go to zero in the steady state. Consequently, the contribution from 
(j • b) mu st b e as large as that of {J ■ B), and of opposite sign, so that the two cancel. This, together 
with eq. (168), allows us immediately to write down an expression for the equilibrium strength of the 
mean field; 

(B')=Mo|(j-b)|/fci, (171) 

which is now valid for both signs of the helicity of the forcing. The residual helicity (Pouquet, Frisch, 
Leorat 1976), 

H,,, = (lj-u) -{j ■b)/po, (172) 

is small in the nonlinear saturated state and nearly vanishing. [We mention that this is also the case 
in the models of Brandenburg & Subramanian (2000).] Furthermore, for forced turbulence with a well 
defi ned f orcing wavenumber the kinetic helicity may be estimated as {uj ■ u) w k[{u'^). Together with 
eq. (171) we have 

(173) 



where B^^ = iiq{pu'^) and _Bcq is the equipartition field strength. Thus, the mean field can exceed (!) the 
equipartition field by the factor {kj /kiY^'^ . This estimate agrees well with the results of the simulations; 
see Brandenburg (2001b). 



Using eqs. (169) and (|170|) together with eqs. (167) and (168) we can rewrite eq. (165) in the form 



-2fjkl{B )+ 277fciMo|(j-b)|, 



(174) 



where we have taken into account the contribution of the small scale current helicity which is of similar 
magnitude as the large scale current helicity. For the magnetic helicity, on the other hand, the small scale 
contribution is negligible, because 



{a-b)\^Po\{3-b)\/k^^f,o\{J-B)\/kf^\{A-B)\ f'j^] «\{A-B 



(175) 



After the saturation at small and intermediate scales the small scale current helicity is approximately 
constant and can be estimated as 



The solution of (174) is given by 



{B ) = eoB, 



1 



(176) 



(177) 



where eo is a coefficient which, in the present model with a well-defined forcing wavenumber, can be 
approximated by ep ~ fcf/fci. 

This is indeed also the limiting behavior found for a^-dynamos with simultaneous a and rj quenching 
of the form 

- 



l + asB /B2q 



l + r,BB/B^^^ 



(178) 



where as = Vb is assumed. Assuming that the magnetic energy density of the mean field, B^ , is 
approximately uniform (which is well satisfied in the simulations) we can obtain the solution B = B{t) 
of eq. (301) in the form 



B'/{l-B'/Bl^Y+^'^^' =Bf^,e 



2\t 



(179) 
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Figure 30: Evoluti on of 
equations using eq. (|179|) 



{B ) for Runs 1-3, 5, and 6, compared with the solution ( |l79|) of the dynamo 



where 



(180) 



is the kinematic growth rate of the dynamo, -Bini is the initial field strength, and -Bgn is the final field 
strength of the large scale field, which is related to as and ris via 



as = Vb 



A 



Bcq 



(181) 



The full derivation is given in Appendix The significance of this result is that it provides an excellent 
fit to the numerical simulations; see figur e |3C| where we present the evolution of {B ) for the different runs 
of Brandenburg (2001b). Equation (179) can therefore be used to extrapolate to astrophysical conditions. 
The time it takes to convert the small scale field generated by the small scale dynamo to a large scale 
field, Toq, increases linear with the magnetic Reynolds number, R^^. Apart from some coefficients of order 
unity the ratio of T^q to the turnover time is therefore just Ry^- For the sun this ratio would be 10^ — 10^". 
However, before interpreting this result further one really has to know whether or not the presence of 
open boundary conditions could alleviate the issue of very long timescales for the mean magnetic field. 
Furthermore, it is not clear whether the long timescales discussed above have any bearing on the cycle 
period in the case of oscillatory solutions. The reason this is not so clear is because for a cyclic dynamo 
the magnetic helicity in each hemisphere stays always of the same sign and is only slightly modulated. It 
is likely that this modulation pattern is advected precisely with the meridional circulation, in which case 
the helicity could be nearly perfectly conserved in a lagrangian frame. This could provide an interesting 
clue for why the solar dynamo is migrating. The relation between meridional circulation and dynamo 
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wave propagation has been advocated by Durney (1995) and Choiidhuri, Schussler, & Dikpati (1995), 
but helicity conservation would strongly lock the two aspects. 

It is clear that virtually all astrophysical bodies are open, allowing for constant loss of magnetic 
helicity. In the case of the sun significant amounts of magnetic helicity are indeed observed at the solar 
surface (Berger & Ruzmaikin 2000). Significant losses of magnetic helicity are particularly obvious in 
the case of accretion discs which are almost always accompanied by strong outfiows that can sometimes 
be collimated into jets. Thus, dynamo action from accretion disc turbulence would be a good candidate 
for clarifying the significance of open boundaries on the nature of the dynamo. Another reason why 
accretion disc turbulence is a fruitful topic for understanding dynamo action is because the shear is 
extremely strong. In the following we discuss some recent progress that has been made in this field. 



7.7 Joule dissipation from mean and fluctuating fields 

In an MHD fiow the mean magnetic Joule dissipation per unit volume is given by 

Qjouic = Wo(J'). (182) 

Whilst in may astrophysical fiows ry may be very small, \J\ can be large so that Qjouie remains finite even 
in the limit ry ^ 0. One example where this is very important is accretion discs, where Joule dissipation 
(together with viscous dissipation) are important in heating the disc. These viscous and resistive processes 
are indeed the only significant sources of energy supply in discs, and yet the luminosities of discs that 
result from the conversion of magnetic and kinetic energies into heat and radiation can be enormous. 
Much of the work on discs involves mean-field modelling, so it would be interesting to see how the Joule 
dissipation, Q predicted from a mean-field model, 

Q,tac = ^tMo(j\ (183) 

relates to the actual Joule dissipation. In figure ^ we show the evolution of actual and mean-field Joule 
dissipation and compare with an estimate for the rate of total energy dissipation, B^^/t, where r is the 
turnover time. Here we have taken into account that rji is 'catastrophically' quenched using the formulae 
of Brandenburg (2001) with the parameters for Run 3. 

There is no reason a priori that the magnetic energy dissipations from the mean-field model should 
agree with the actual one. It turns out that the mean-field dissipation is a fourth of the actual one, so it 
is definitely significant. It would therefore be interesting so see how those two dissipations compare with 
each other in other models. 



7.8 Possible pitfalls in connection with hyperresistivity 

In many astrophysical applications hyperresistivity and hyperviscosity are sometimes used in order to 
concentrate the effects of magnetic diffusion and viscosity to the smallest possible scale. The purpose of 
this section is to highlight possible spurious artifacts associated with this procedure. As we have seen 
above, large scale dynamos can depend on the microscopic magnetic diffusivity and must therefore be 
affected when it is replaced by hyperresistivity. The resulting modifications that are to be expected are 



easily understood: on the right hand side of eq. (165) the term (J-B) needs to be replaced by ((V A)-B 



This leads to a change of the relative importance of small and large scale contributions, which therefore 
changes eq. (FtsI) to 



t2, 



(B'>=(^^j <• (184) 

Thus, the final field strength will be even larger than before: instead of a factor of 5 superequipartition 
(for fcf = 5) one now expects a factor of 125. Recent simulations by Brandenburg & Sarson (2001) have 
indeed confirmed this tendency. The main conclusion is that hyperresistivity can therefore be used to 
address certain issues regarding large magnetic Reynolds numbers that are otherwise still inaccessible. 
On the other hand, the results are in some ways distorted and need therefore be interpreted carefully. 
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Figure 31: Joule dissipation for Run 3 (solid line), compared with the Joule dissipation estimated for 
a corresponding mean-field model (dashed line). An estimate for the rate of total energy dissipation, 
-Bg /r, is also given. 



7.9 Remarks on accretion disc turbulence 

We have already mentioned the possibility of dynamo action in accretion discs. Accretion discs have been 
postulated some 30 years ago in order to explain the incredibly high luminosities of quasars. Only in the 
past few years has direct imaging of accretion discs become possible, mostly due to the Hubble Space 
Telescope. Accretion discs form in virtually all collapse processes, such as galaxy and star formation. In 
the latter case the central mass is of the order of one solar mass, while in the former it is around 10^ solar 
masses and is concentrated in such a small volume that that it must be a black hole. If the surrounding 
matter was nonrotating, it would fall radially towards the center. But this is unrealistic and even the 
slightest rotation relative to the central object would become important eventually as matter falls closer 
to the center. 

If there was no effective diffusive process in discs, the angular momentum of the matter would stay 
with the gas parcels, and since the gravity force is balanced by the corresponding centrifugal force, the 
gas would never accrete. However, the angular velocity of the gas follows a r~^^^ Kepler law, so the 
gas is differentially rotating and one may expect shear instabilities to occur that would drive turbulence 
and hence turbulent dissipation. Unfortunately, however, the story is not so simple. Discs are both 
linearly stable (Stewart 1975) and probably also nonlinearly stable (Hawley, Gammie, & Balbus 1996). 
Nevertheless, in the presence of a magnetic field there is a powerful linear instability (Balbus & Hawley 
1991), and subsequent work has shown that this instability is indeed capable of driving the instability 
and hence turbulence. 

One of the key outcomes of such simulations is the rate of turbulent dissipation, which determines the 

rate of angular momentum transport and correspondingly the rate at which orbital kinetic and potential 
energy is liberated in the form of heat. This is normally expressed in terms of a turbulent viscosity (e.g. 
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Figure 32: Dependence of ass on time and the mean magnetic field strength for a local accretion disc 
model (Run B of Brandenburg et al. 1996a). Here i?o = (MopCg)^/^ is the thermal equipartition field 
strength and Trot — 2ti/VI the local rotation period. In the left hand panel the dotted line represents 
the actual data and the solid line gives the fit obtained by correlating ass with the mean magnetic field 
(right hand panel). 



Frank, King, & Raine 1992), but it may equally well be expressed in terms of the horizontal components 
of the Reynolds and Maxwell stress tensors. The stress may then be normalized by {pc^) to give a 
nondimensional measure (called ass) of the ability of the turbulence to transport angular momentum 
outward (if ass > 0)- This ass is indeed always positive, see figure |3^, but it fluctuates significantly about 
a certain mean value. These fluctuations are in fact correlated with the energy in the mean magnetic 
field, (-B)^, as is shown in the right hand panel of figure ^ This mean magnetic field shows regular 
reversals combined with a migration away from the midplane, as can be seen in figure 

The evolution of the mean magnetic field found in the simulations is reminiscent of the behavior 
known from mean-field a — dynamos. Further details regarding this correspondence (relation between 
the value of a and cycle period, field parity for different boundary conditions, etc.) can be found in recent 
reviews of the subject (e.g., Brandenburg 1998, 2000). 



7.10 Connection with the solar dynamo problem 

The disc simulations have shown that a global large scale field can be obtained even in cartesian geometry. 
The detailed behavior of this large scale field depends of course on the boundary conditions adopted 
(Brandenburg 1998), and will therefore be different in different geometries. Nevertheless, the very fact 
that large scale dynamo action is possible already in simple cartesian geometry is interesting. 

In figure |^ we show the evolution of magnetic and kinetic energies as well as the magnitudes of 
the large scale field for a simulation of a convectively driven dynamo in the presence of large scale shear 
(Brandenburg, Nordlund, & Stein 2001). It turns out that the ratio of the magnetic energies in large scale 
fields relative to the total field, (-B)^/ (B^), which is a measure of the filling factor of the magnetic field, is 
around 15% when the field has reached saturation, i.e. when the field growth has stopped. This is similar 
to the case of isotropic nonmirror-symmetric turbulence considered in § On the average, however, 
the magnetic field is then directed into the negative j/-direction (corresponding to the negative azimuthal 
direction in spherical geometry) , but there is a weak and more noisy cross-stream field component directed 
in the positive x-direction (pointing north). 

When the field has reached saturation, the mean field direction is approximately constant. Although 
the magnitude of this mean field fluctuates somewhat, the sign is always the same. Thus, this simulation 
shows no cycles, which are so characteristic of the solar dynamo. However, since those features, including 
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Figure 33: Butterfly (space- time) diagram of the poloidal and toroidal magnetic field components 

averaged over the two horizontal (x and y) directions for the local accretion disc model of (Brandenburg 
et al. 1996a, Run B). Note that the poloidal field is much more noisy than the toroidal field, and that 
there is a clear outward migration of magnetic field. 

the field geometry depend strongly on boimdary condition and on the location of the boundary conditions, 
this disagreement is to be expected, and one would really need to resort to global simulations in spherical 
geometry. 
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Figure 34: Evolution of several quantities for a convective dynamo model with shear: kinetic and 
magnetic energies (dotted and solid lines; first panel), mean latitudinal and toroidal fields (dotted and 
solid lines; second panel), mean magnetic field in a linear scale (third panel), and the filling factor (fourth 
panel). Energies and mean magnetic fields are given in units of the equipartition value, B^q = 47r(M^). 



7.11 Dynamics of the overshoot layer 

Late-type stars with outer convection zones have an interface between the convection zone proper and 
the radiative interior. This leads to some additional dynamics that is important to include, especially in 
connection with the dynamo problem. This interface is the layer where magnetic flux can accumulate, 
i.e. not necessarily the layer where the dynamo operates; see the discussion in Brandenburg (1994). The 
accumulation is a consequence of turbulent pumping down the turbulence intensity gradient and the effect 
was seen clearly in video animations reported by Brandenburg & Tuominen (1991) and was analyzed in 
detail by Nordlund et al. (1992). Tobias et al. (1998) have studied the effect in isolation starting with an 
initial magnetic field distribution as opposed to a dynamo-generated field. 

The flow dynamics changes drastically as one enters this overshoot layer. The stabilizing buoyancy 
effect provides a restoring force on a downward moving element, which can give rise to gravity waves 
that could be driven by individual plumes. This leads to a marked wavy pattern that can extend deep 
into the lower overshoot layer, as seen in flgure ^ where we have plotted the vertical rms velocity as a 
function of depth and time. These waves extend a major fraction into the stably stratified layer beneath 
the convection zone, but are damped eventually. The typical period of such events is seen to be around 
20 (in units of d/g), where d is the depth of the unstable (convective) layer. This is comparable with 
the mean Brunt- Vaisala frequency, 

^lv--9-V(,s/cp). (185) 
which is around 0.3 in the overshoot layer; see figure ^ 
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Figure 35: Space-time diagram of the vertical rms velocities for a nonmagnetic convection zone model of 
Brandenburg et al. (2001). Note the propagation of isolated plumes in more or less regular time intervals. 
Note also that the wavy pattern extends well into the convection zone proper (0.5 < z < 1), and that 
the plumes appear to propagate at an approximately constant speed towards the bottom. This speed is 
around 0.1, which is comparable to the rms velocity in the runs. 

Technically the presence of a lower overshoot layer provides a formidable challenge, because there the 
dynamics is governed by the very slow thermal timescale. This can lead to problems if the properties of 
the upper surface layers change which could affect the entropy in deeper parts of the convection zone. 
This in turn will also affect the stratification in the rest of the radiative interior. Since this can only 
happen on a thermal timescale (there is no turbulence in these layers to speed things up) it takes very 
long before one arrives at a new statistically stable state. This problem is of course also encountered when 
starting with initial conditions that are derived under too unrealistic assumptions. There are essentially 
two different approaches to this: one either considers a toy model where dynamical and thermal timescales 
are artificially brought closer together or one adopts an implicit code which allows the use of somewhat 
longer timesteps. The former approach implies that one adopts input fluxes that exceed those in real 
stars, but the good news is that turbulent velocities and temperature fluctuations vary with changing 
flux in exactly the way that is expected from mixing length theory; see Brandenburg, Nordlund, & Stein 
(2001) for details. 

8 Conclusions 

Many phenomena in astrophysics show direct manifestations of turbulence. As in the case of accretion 
discs, without turbulence there would be no enhanced dissipation, no heating of the disc, and hence no 
emission. Magnetic fields are of major importance as is evidenced again by the example of accretion 
discs, where the turbulence is a direct consequence of the presence of magnetic fields (Balbus-Hawley 
instability). Although magnetic fields in discs could have primordial origin and could just have been 
compressed during their formation, it is also clear that discs are actually a favorable candidate for 
producing strong large scale magnetic fields, as shown by the local simulations discussed above. 
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Figure 36: Modulus of the Brunt- Vaisala frequency for a run with polytropic index m = —0.9 and a 
nondimensional input flux JFconv = 0.01. The various curves are for different values of the radiative flux, 
but fixed nominal convective flux. Small values of J-"rad/.^tot are typical for the upper layers of the solar 
convection zone. 

Other bodies where strong dynamo action is possible arc stars. Simulations of convection have shown 
that strong small scale magnetic fields are possible (Meneguzzi & Pouquet 1989, Nordlund et al. 1992, 
Brandenburg et al. 1996b, Cattaneo 1999), but if there is strong shear an intense large scale field can 
also emerge. 

The use of nonconservative high order schemes has proved useful in all those applications. They 
are easy to implement and to modify, but they are also reasonably accurate. In this chapter we have 
illustrated the behavior of such schemes using various test problems. Using potential enthalpy and entropy 
as the main thermodynamic variables has a number of advantages, especially in connection with strongly 
stratified flows near a central object with a deep potential well, which is relevant to studying outflow 
phenomena. Contrary to common belief, high order schemes arc not more vulnerable to Gibbs phenomena 
near discontinuities. Instead, in simple advection tests high order methods are able to produce smoother 
solutions with less viscosity, which is important for accurate modeling of turbulence. 

In the last part of this chapter we have briefly mentioned some astrophysical applications of simula- 
tions using high order schemes where hydromagnetic turbulence played an important role. In the next 
few years we may expect a dramatic increase in the quality and predictive power of such simulations, 
as larger computers become available. Already now a number of very promising simulations arc emerg- 
ing. There is important work addressing the stability of astrophysical jets in three dimensions (Ouyed, 
Clarke, & Pudritz 2000). Also worthwhile mentioning are recent high resolution simulations by Hawley 
(2000) of three-dimensional accretion tori in global geometry. What remains to be done in this fleld 
is a proper connection between disc physics and the launching mechanism of jets. This would require 
incorporating proper thermodynamics allowing for radiative cooling and magnetic heating in particular. 
Global simulations would also be highly desirable to address the global stellar dynamo problem. It would 
be interesting to see how the dynamo works in fully convective stars, for example. This problem is in 
some ways simpler than the solar dynamo problem, because one does not need to worry about the lower 
overshoot layer where the relevant timescales are much longer than in the convection zone proper. 

Acknowledgements: I am indebted to Ake Nordlund for teaching me many of the methods and 




52 



techniques that are reflected to some extend in the present work. I thank Wolfgang Dobler and Petri 
Kapyla for their careful reading of the manuscript and for pointing out a number of mistakes in the 
manuscript. 

References 

Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001 A staggered grid, Lagrangian- 

Eulerian remap code for 3D MHD simulations. J. Comp. Phys. 171, 151-181. 
Balbus, S. A. & Hawley, J. F. 1991 A powerful local shear instability in weakly magnetized disks. I. 

Linear analysis. Astrophys. J. 376, 214-222. 
Balsara, D., & Pouquet, A. 1999 The formation of large-scale structures in supersonic magnetohydro- 

dynamic flows. Phys. Plasmas 6, 89-99. 
Benson, D. J. 1992 Computational methods in Lagrangian and Eulerian hydrocodes. Comp. Meth. 

Appl. Mech. Eng. 99, 235 394. 
Berger, M. A., & Ruzmaikin, A. 2000 Rate of helicity production by solar rotation. J. Geophys. Res. 

105, 10481-10490. 

Biskamp, D. 1994 Cascade models for magnetohydro dynamic turbulence. Phys. Rev. E50, 2702 2711. 
Biskamp, D., & Miillcr, W.-C. 1999 Decay laws for three-dimensional magnetohydrodynamic turbulence. 

Phys. Rev. Lett. 83, 2195 2198. 
Bigazzi, A. 1999 Models of small-scale and large-scale dynamo action. PhD thesis, Universita Degli 

Studi Dell'Aquila. 

Blackman, E. C, & Field, G. F. 2000 Constraints on the magnitude of a in dynamo theory. Astrophys. 
J. 534, 984-988. 

Blandford, R. D. & Payne, D. G. 1982 Hydromagnetic flows from accretion discs and the production of 

radio jets. Monthly Notices Roy. Astron. Soc. 199, 883-903. 
Brandenburg, A. 1992 Energy spectra in a model for convective turbulence. Phys. Rev. Lett. 69, 

605-608. 

Brandenburg, A. 1994 Solar Dynamos: Computational Background. In Lectures on Solar and Planetary 

Dynamos (ed. M. R. E. Proctor & A. D. Gilbert), pp. 117-159. Cambridge University Press. 
Brandenburg, A. 1998 Disc Turbulence and Viscosity. In Theory of Black Hole Accretion Discs (ed. M. 

A. Abramowicz, G. Bjornsson & J. E. Pringlc), pp. 61 86. Cambridge University Press. 
Brandenburg, A. 1999 Helicity in large-scale dynamo simulations. In Magnetic Helicity in Space and 

Laboratory Plasmas (ed. M. R. Brown, R. C. Canfield, A. A. Pevtsov), pp. 65-73. Geophys. 

Monograph 111, American Geophysical Union, Florida. 
Brandenburg, A. 2000 Dynamo-generated turbulence and outflows from accretion discs. Phil. Trans. 

Roy. Soc. Lond. A 358, 759-776. 
Brandenburg, A. 2001a Magnetic mysteries. Science 292, 2440-2441. 

Brandenburg, A. 2001b The inverse cascade and nonlinear alpha-effect in simulations of isotropic helical 

hydromagnetic turbulence. Astrophys. J. 550, 824-840. 
Brandenburg, A., & Donner, K. J. 1997 The dependence of the dynamo alpha on vorticity. Monthly 

Notices Roy. Astron. Soc. 288, L29-L33. 
Brandenburg, A., & Schmitt, D. 1998 Simulations of an alpha-effect due to magnetic buoyancy. Astron. 

Astrophys. 338, L55-L58. 

Brandenburg, A., & Subramanian, K. 2000 Large scale dynamos with ambipolar diffusion nonlinearity. 
Astron. Astrophys. 361, L33-L36. 

Brandenburg, A., Tuominen, I. 1991 The solar dynamo. In The Sun and cool stars: activity, magnetism, 
dynamos, lAU Coll. 130 (ed. I. Tuominen, D. Moss & G. Riidiger), pp. 223-233. Lecture Notes 
in Physics 380, Springer- Verlag. 

Brandenburg, A., Enqvist, K., Olesen, P. 1996 Large-scale magnetic fields from hydromagnetic turbu- 
lence in the very early universe. Phys. Rev. D 54, 1291-1300. 



53 



Brandenburg, A., Nordlund, A., & Stein, R. F. 2000 Astrophysical convection and dynamos. In 
Geophysical and Astrophysical Convection (ed. P. A. Fox & R. M. Kerr), pp. 85-105. Gordon and 
Breach Science Publishers. 

Brandenburg, A., Nordlund, A., & Stein, R. F. 2001 "Simulation of a convective dynamo with imposed 
shear," Astron. Astrophys. (to be submitted). 

Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995 Dynamo generated turbulence and 
large scale magnetic fields in a Keplerian shear flow. Astrophys. J. 446, 741-754. 

Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1996a The disk accretion rate for dynamo 
generated turbulence. Astrophys. J. Letters 458, L45-L48. 

Brandenburg, A., Jennings, R. L., Nordlund, A., Rieutord, M., Stein, R. F., & Tuominen, I. 1996b Mag- 
netic structures in a dynamo simulation. J. Fluid Mech. 306, 325-352. 

Brandenburg, A., Dobler, W., Shukurov, A., & von Rekowski, B. 2000 "Pressure-driven outflow from a 
dynamo active disc," Astron. Astrophys. (submitted). astro-ph/0003174| 

Brandenburg, A., Chan, K. L., Nordlund, A., Stein, R. F. 2001 "The effect of the radiative background 
flux in convection simulations," Astron. Astrophys. (to be submitted). 

Brand enburg A., fc Sarson, G. R. 2001 The effect of hyperdiffusivity on turbulent dynamos with helicity. 



http : 1 1 www ■ nordit a ■ dk/ ~br eaidenb/tmp/ graeme/paper . ps . g; 



Campbell, C. G. 1999 Launching of accretion disc winds along dynamo-generated magnetic fields. 

Monthly Notices Roy. Astron. Soc. 310, 1175-1184. 
Campbell, C. G. 2000 An accretion disc model with a magnetic wind and turbulent viscosity. Monthly 

Notices Roy. Astron. Soc. 317, 501-527. 
Canuto, C, Hussaini, M. Y., Quarteroni, A., & Zang, T. A. 1988 Spectral Methods in Fluid Dynamics. 

Springer, Berlin. 

Cattaneo, F., 1999 On the origin of magnetic fields in the quiet photosphere. Astrophys. J. 515, 
L39-L42. 

Cattaneo, F., & Hughes, D. W. 1996 Nonlinear saturation of the turbulent alpha effect. Phys. Rev. E 
54, R4532-R4535. 

Chan, K. L., & Sofia, S. 1986 Turbulent compressible convection in a deep atmosphere. II. Tests on the 
validity and limitation of the numerical approach. Astrophys. J. 307, 222-241. 

Chan, K. L., & Sofia, S. 1989 Turbulent compressible convection in a deep atmosphere. III. Results of 
three-dimensional computations. Astrophys. J. 336, 1022-1040. 

Choudhuri, A. R., Schiissler, M., & Dikpati, M. 1995 The solar dynamo with meridional circulation. 
Astron. Astrophys. 303, L29-L32. 

Collatz, L. 1966 The numerical treatment of differential equations. Springer- Verlag, New York., p. 164. 

Colella, P., & Woodward, P. R. 1984 The piecewise parabolic method (PPM) for gas-dynamical simu- 
lations. J. Comp. Phys. 54, 174-201. 

Dobler, W., Brandenburg, A., & Shukurov, A. 1999 Pressure-driven outflow and magneto-centrifugal 
wind from a dynamo active disc. In Plasma Turbulence and Energetic Particles in Astrophysics 
(ed. M. Ostrowski & R. Schlickeiser) , pp. 347-352. Publ. Astron. Obs. Jagiellonian Univ., Cracow. 

Dodd, R. K., Eilbeck, J. C, Gibbon, J. D., & Morris, H. C. 1982 Solitons and nonlinear wave equations. 
Academic Press: London. 

Durney, B. R. 1995 On a Babcock-Leighton dynamo model with a deep-seated generating layer for the 

toroidal magnetic fleld. II. Solar Phys. 166, 231-260. 
Fox, P. A., Theobald, M. L., & Sofia, S. 1991 Compressible magnetic convection: formulation and 

two-dimensional models. Astrophys. J. 383, 860-881. 
Frank, J., King, A. R., & Raine, D. J. 1992 Accretion power in astrophysics. Cambridge: Cambridge 

Univ. Press. 

Frisch, U., Pouquet, A., Leorat, J., Mazure, A. 1975 Possibility of an inverse cascade of magnetic helicity 

in hydrodynamic turbulence. J. Fluid Mech. 68, 769-778. 
Galsgaard, K. & Nordlund, A. 1996 Heating and activity of the solar corona: I. boundary shearing of 



54 



an initially homogeneous magnetic-field. /. Geophys. Res. 101, 13445-13460. 
Galsgaard, K. & Nordlund, A. 1997a Heating and activity of the solar corona: II. Kink instability in a 

flux tube. J. Geophys. Res. 102, 219-230. 
Galsgaard, K. & Nordlund, A. 1997b Heating and activity of the solar corona: HI. Dynamics of a low 

beta plasma with three-dimensional null points. /. Geophys. Res. 102, 231-248. 
Glatzmaier, G. A., & Roberts, P. H. 1995 A three-dimensional self-consistent computer simulation of a 

geomagnetic field reversal. Nature 377, 203-209. 
Glatzmaier, G. A., & Roberts, P. H. 1996 Rotation and magnetism of Earth's inner core. Science 274, 

1887-1891. 

Grauer, R., Marliani, C., & Germaschewski, K. 1998 Adaptive mesh refinement for singular solutions 
of the incompressible Euler equations. Phys. Rev. Lett. 80, 4177-4180. 

Greenhill, L. J., Gwinn, C. R., Schwartz, C., Moran, J. M., & Diamond, P. J. 1998 Coexisting conical 
bipolar and equatorial outflows from a high-mass protostar. Nature 396, 650-653. 

Hawley, J. F. 2000 Global magnetohydrodynamical simulations of accretion tori. Astrophys. J. 528, 
462-479. 

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995 Local three-dimensional magnetohydrodynamic 

simulations of accretion discs. Astrophys. J. 440, 742-763. 
Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996 Local three dimensional simulations of an accretion 

disk hydromagnetic dynamo. Astrophys. J. 464, 690-703. 
Hodapp, K.-W. & Ladd, E. F. 1995 Bipolar jets from extremely young stars observed in molecular 

hydrogen emission. Astrophys. J. 453, 715-720. 
Ji, H., Prager, S. C., Almagri, A. F., Sarff, J. S., & Toyama, H. 1996 Measurement of the dynamo effect 

in a plasma. Phys. Plasmas 3, 1935-1942. 
Kerr, R. M., & Brandenburg, A. 1999 Evidence for a singularity in ideal magnetohydrodynamics: 

implications for fast reconnection. Phys. Rev. Lett. 83, 1155-1158. 
Kleeorin, N. I, Moss, D., Rogachevskii, I., & Sokoloff, D. 2000 Helicity balance and steady-state strength 

of the dynamo generated galactic magnetic fleld. Astron. Astrophys. 361, L5-L8. 
Korpi, M. J., Brandenburg, A., Shukurov, A., Tuominen, I., & Nordlund, A. 1999 A supernova regulated 

interstellar medium: simulations of the turbulent multiphase medium. Astrophys. J. Letters 514, 

L99-L102. 

Lele, S. K. 1992 Compact flnite difference schemes with spectral-like resolution. J. Comp. Phys. 103, 
16-42. 

LeVeque, R. J., Mihalas, D., Dorfl, E. A., & Miiller, E. 1998 Computational methods for astrophysical 

fluid flow. Springer, Berlin. 
Lohse, D., & Miiller-Groeling, A. 1995 Bottleneck effects in turbulence: Scaling phenomena in r versus 

p space. Phys. Rev. Lett. 74, 1747-1750. 
Mac Low, M.-M., Klessen, R. S., & Burkert, A. 1998 Kinetic energy decay rates of supersonic and 

super-Alfvenic turbulence in star-forming clouds. Phys. Rev. Lett. 80, 2754-2757. 
Meneguzzi, M., & Pouquet, A. 1989 Turbulent dynamos driven by convection. J. Fluid Mech. 205, 

297-312. 

Nordlund, A. 1982 Numerical Simulations of the Solar Granulation I. Basic Equations and Methods. 

Astron. Astrophys. 107, 1-10. 
Nordlund, A. 1985 Solar convection. Solar Phys. 100, 209-235. 



Nordlund, A., fc Galsgaard, K. 1995 A 3D MHD code for Parallel Computers, http : / / www . astro . ku 



dk/~aake/NumericalAstro/papers/kg/mlid . ps . gz . 



Nordlund, A., & Stein, R. F. 1990 3-D Simulations of Solar and Stellar Convection and Magnetocon- 

vection. Comput. Phys. Commun. 59, 119-125. 
Nordlund, A., Galsgaard, K., & Stein, R. F. 1994 Magnetoconvection and Magnetoturbulence. In Solar 

surface magnetic fields (ed. R. J. Rutten & C. J. Schrijver), pp. 471-498. NATO ASI Series, Vol. 

433. 



55 



Nordlund, A., Brandenburg, A., Jennings, R. L., Rieutord, M., Ruokolainen, J., Stein, R. F., & Tnominen, 
I. 1992 Dynamo action in stratified convection with overshoot. Astrophys. J. 392, 647 652. 

Ouyed, R., Pudritz, R. E. & Stone, J. M. 1997 Episodic jets from black holes and protostars. Nature 
385, 409-414. 

Ouyed, R., & Pudritz, R. E. 1997a Numerical simulations of astrophysical jets from keplerian discs. I. 

Stationary models. Astrophys. J. 482, 712-732. 
Ouyed, R., & Pudritz, R. E. 1997b Numerical simulations of astrophysical jets from keplerian discs. II. 

Episodic outflows. Astrophys. J. 484, 794-809. 
Ouyed, R., & Pudritz, R. E. 1999 Numerical simulations of astrophysical jets from keplerian discs. III. 

The effects of mass loading. Monthly Notices Roy. Astron. Sac. 309, 233-244. 
Ouyed, R., Clarke, D. A. & Pudritz, R. E.: 2000, "3-Dimensional simulations of astrophysical jets from 

keplerian accretion Disks I: stability issues," Astrophys. J. (in press, scheduled for the 
issue)Padoan, P., Nordlund, A., & Jones, B. J. T. 1997 The universality of the stellar mass function. 

Astrophys. J. 288, 145 152. 
Padoan, P., & Nordlund, A. 1999 A super-Alfvnic model of dark clouds. Astrophys. J. 526, 279-294. 
Passot, T., & Pouquet, A. 1987 Numerical simulation of compressible homogeneous flows in the turbulent 

regime. J. Fluid Mech. 181, 441-466. 
Pelletier, G., & Pudritz, R. E. 1992 Hydromagnetic disk winds in young stellar objects and active 

galactic nuclei. Astrophys. J. 394, 117-138. 
Peterkin, R. E., Frese, M. H., & Sovinec, C. R. 1998 Transport of magnetic flux in an arbitrary coordinate 

ALE code. J. Comp. Phys. 140, 148-171. 
Porter, D. H., Pouquet A., & Woodward, P. R. 1992 Three-dimensional supersonic homogeneous 

tiirbulenee: a numerical study. Phys. Rev. Lett. 68, 3156-3159. 
Porter, D. H., Pouquet A., & Woodward, P. R. 1994 Kolmogorov-like spectra in decaying 2-dimensional 

supersonic flows. Phys. Fluids 6, 2133-2142. 
Pouquet, A.. Frisch, U., Leorat, J. 1976 Strong MHD helical turbulence and the nonlinear dynamo 

effect. J. Fluid Mech. 77, 321 354. 
Rast, M. P., Nordlund, A., Stein, R.F., & Toomre, J. 1993 Ionization effects in three-dimensional 

granulation simulations. Astrophys. J. Letters 408, L53-L56. 
Rast, M. P., & Toomre, J. 1993a Compressible convection with ionization. I. Stability, flow asymmetries, 

and energy transport. Astrophys. J. 419, 224-239. 
Rast, M. P., & Toomre, J. 1993b Compressible convection with ionization. II. Thermal boundary-layer 

instability. Astrophys. J. 419, 240-254. 
Rekowski, M. v., Riidiger, G., & Elstner, D. 2000 Structure and magnetic configurations of accretion 

disk-dynamo models. Astron. Astrophys. 353, 813-822. 
Roettiger, K., Stone, J. M., & Burns, J. O. 1999 Magnetic field evolution in merging clusters of galaxies. 

Astrophys. J. 518, 594-602. 
Rognvaldsson, O. E., Nordlund, A., & Sommer-Larsen, J. 2001 Cooling flows and disk galaxy formation. 

preprint. 

Rogallo, R. S. 1981 Numerical experiments in homogeneous turbulence. NASA Tech. Memo. 81315. 
Sanchez- Salcedo, F. J., & Brandenburg, A. 1999 Deceleration by dynamical friction in a gaseous medium. 

Astrophys. J. Letters 522, L35-L38. 
Sanchez-Salcedo, F. J., & Brandenburg, A. 2001 Dynamical friction of bodies orbiting in a gaseous 

sphere. Monthly Notices Roy. Astron. Soc. 322, 67-78. 
Sod, G. A. 1978 A survey of several finite difference methods for systems of nonlinear hyperbolic 

conservation laws. J. Comp. Phys. 27, 1 31. 
Stanescu, D. & Habashi, W. G. 1998 2A^-storage low dissipation and dispersion Runge-Kutta schemes 

for computational acoustics. J. Comp. Phys. 143, 674-681. 
Steffen, M., Ludwig, H.-G., & KriiB, A. 1989 A numerical study of solar granular convection in cells of 

different horizontal dimensions. Astron. Astrophys. 123, 371-382. 



56 



Stein, R.F., & Nordlund, A. 1989 Topology of convection beneath the solar surface. Astrophys. J. 
Letters 342, L95-L98. 

Stein, R.F., & Nordlund, A. 1998 Simulations of solar granulation. I. General properties. Astrophys. J. 
499, 914-933. 

Stone, J. M., Norman, M. 1992a ZEUS-2D: A radiation magnetohydrodynamics code for astrophysical 
flows in two space dimensions; I. The hydrodynamic algorithms and tests. Astrophys. J. 80, 
753-790 

Stone, J. M., Norman, M. 1992b ZEUS-2D: A radiation magnetohydrodynamics code for astrophysical 
flows in two space dimensions: II. The magnetohydrodynamic algorithms and tests. Astrophys. J. 
80, 791-818 

Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996 Three dimensional magnetohydrody- 

namical simulations of vertically stratified accretion disks. Astrophys. J. 463, 656-673. 
Stewart, J. M. 1975 The hydrodynamics of accretions discs 1: Stability. Astron. Astrophys. 42,95-101. 
Suzuki, E. & Toh, S. 1995 Entropy cascade and temporal intermittency in a shell model for convective 

turbulence. Phys. Rev. E 51, 5628-5635. 
Thelen, J.-C. 2000 A mean electromotive force induced by magnetic buoyancy instabilities. Monthly 

Notices Roy. Astron. Soc. 315, 155-164. 
Tobias, S. M., Brummell, N. H., Clune, T. L., & Toomre, J. 1998 Pumping of magnetic fields by 

turbulent penetrative convection. Astrophys. J. Letters 502, L177-L177. 
Toh, S., & lima, M. 2000 Dynamical aspect of entropy transfer in free convection turbulence. Phys. 

Rev. E 61, 2626-2639. 

Urpin, v., & Brandenburg, A. 1998 Magnetic and vertical shear instabilities in accretion discs. Monthly 

Notices Roy. Astron. Soc. 294, 399-406. 
Vainshtein, S. I., & Cattaneo, F. 1992 Nonlinear restrictions on dynamo action. Astrophys. J. 393, 

165-171. 

Vishniac, E. T., & Brandenburg, A. 1997 An incoherent a — f2 dynamo in accretion disks. Astrophys. 
.J. 475, 263-274. 

Volk, H. J. & Atoyan, A. M. 1999 Clusters of galaxies: magnetic fields and nonthermal emission. 

Astroparticle Phys. 11, 73-82. 
Williamson, J. H. 1980 Low-storage Runge-Kutta schemes. J. Comp. Phys. 35, 48-. 
Ziegler, U., & Riidiger, G. 2000 Angular momentum transport and dynamo-effect in stratified, weakly 

magnetic disks. Astron. Astrophys. 356, 1141-1148. 



Appendix 

A Centered, onesided and semi-onesided derivatives 

In §^ we gave the centered finite difference formulae for schemes of different order. Here we first describe 
the method for determining the finite difference formulae for second order, but the generalization to higher 
order is straightforward. We also give the corresponding expressions for one-sided and semi-onesided finite 
difference formulae. 

We want to write the derivative f'{x) at the point 

/i = a-i/j-i + oo/o + ai./i, (186) 

where // = f'{xi), fi-i — f{xi — 5x)^ and /i+i = f{xi + 5x). To determine the coefficients a_i, ag and 
ai we expand f{x) up to second order 

/(x) = Co + cix + C2a;^. (187) 

The first derivative is then 

f'{x) = ci + 2c2X. (188) 
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In particular, the value at x = is just /'(O) = ci. Likewise, we have /"(O) — 2ci To determine all 
coefficients we make use of our knowledge at the neighboring points around x^, i.e. we use the function 
values /(xi — bx) = /i-i, j[xi) = fi, and f{xi + Sx) = fi+i, so we have 



fi-l 

h 

fi+l 



This can be written in matrix form 




= Cf) + ci{~Sx) + C2{—6x)'^ , 
= Co, 

= co + ci{+Sx) + C2i+Sxy. 



•(-!)" (-1)1 (-1)2 
0° Qi 02 
1° li 12 




(where 



1), or 



and so we obtain the coefficients as 



f = Mc, 



(189) 
(190) 
(191) 



(192) 



(193) 



(194) 



To calculate /' we need the value of ci, see eq. (188), and so the coefficients a„ needed to express the 
derivative are a_i = (7W~i)io, oo = (-^^^)iij ^^nd ai = {M^^)i2, i.e. all points of the inverted matrix 
in the second row. The resulting formula for f- is well-known, 



/; = (-/.-i + /.+i)/(25x). 

The corresponding result for the second derivative is 



(195) 



(196) 



On the boundaries we have to calculate for derivative using only points inside the domain, which is 
explained in the next subsection for second order accuracy, but again the generalization to higher order 
is straightforward and only the results will be listed. 



A.l One-sided 2nd order derivatives 

Again, we want to write the derivative f'{x) as 

/;-ao/o + 01/1+02/2, (197) 

but now 

/. \ /O" Qi 02\ / CO 
/,+i = 1" 1^ 12 cife I (198) 
V2" 21 22/ \c25x^ 

Thus, one arrives at 

/; = (-3/. + 4/,+i - /^+2)/(2<5x). (199) 
Correspondingly, for the second derivative we have 

/r = (2/, - 5/,+i + 4/,+2 - /.+3)/fe2. (200) 
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A. 2 4th order derivatives 

First derivatives 

/; = (/,-2 - 8/,_i + 8/,+i - .h+2)/il2dx) (201) 

- (-3/.-1 - 10/. + 18/,+i - 6f,+2 + /,+3)/(125x) (202) 

/; = (-25/, + 48/.+1 - 36/,+2 + 16/.+3 - 3/.+4)/(12fe) (203) 

Second derivatives 

f- = {-h-2 + 16/.-1 - 30/, + 16/,+i - /,+2)/(12fe2) (204) 

/;' = (ll/,_i - 20/, + 6/,+i + 4/,+2 - h+z)/{l25x'') (205) 

/r = (35/, - 104/,+i + 114/,+2 - 56/,+3 + ll/,+4)/(12(5a;2) (206) 

A. 3 6th order derivatives 

First derivatives 

/; = (-/,-3 + 9/,-2 - 45/,_i + 45/,+i - 9/,+2 + /,+3)/(60fa) (207) 

/; - (2/,_2 - 24/,_i - 35/, + 80/,+i - 3G/,+2 + 8/,+3 - /,+4)/(605x) (208) 

/; - (-10/,_i - 77/, + 150/,+! - 100/,+2 + 50/,+3 - 15/,+4 + 2/,+5)/(60^a;) (209) 

/; = (-147/, + 360/,+! - 450/,+2 + 400/,+3 - 225/,+4 + 72/,+5 - 10/,+6)/(60fe) (210) 
Second derivatives 

f[' = (2/,_3 - 27/,_2 + 27G/,_i - 490/, + 27G/,+i - 27/,+2 + 2/,+3)/(180fa2) (211) 

/f = (-13/,_2 + 228/,_i - 420/, + 200/,+i + 15/,+2 - 12/,+3 + 2/,+4)/(180fe2) (212) 

/," = (137/,_i - 147/, - 255/,+! + 470/,+2 - 285/,+3 + 93/,+4 - 13/,+5)/(180fc2) (213) 

/," = (812/, - 3132/,+! + 5265/,+2 - 5080/,+3 + 2970/,+4 - 972/,+5 + 137/,+6) /{imx") (214) 



B The 2iV-RK3 scheme 

If N is the number of variables to be updated from one timestep to the next, the 2Af-schemes require only 
2 X N variables to be stored in memory at any time. This is better than for the standard Runge-Kutta 
schemes. The general iteration formula is 



(215) 



For a third order scheme we have i = 1, 3. In order to advance the variable u from u^"^ at time t'") to 
at time t^"+^^ = t^"^ + ft. we set in eq. (|l|) 



wo = M^"^ and W3 = 



(216) 



with ui and U2 being intermediate steps. In order to be able to calculate the first step, « = 1, for which 
no Wi-i = wq exists, we hav e to require ai — 0. Thus, we are left with 5 unknowns, a2, as, /3i, /32, and 
P3. We write down eq. (^151) in explicit form for i = 1, ...,3, 



wi ^ hF{to,UQ), ui^UQ+PiWi. 

W2 — a2Wi + hF{ti,Ui), U2 = Ui + (32W2- 
W3 = a3W2 + hF{t2,U2), U3 ^ U2 + (33W3. 



(217) 

(218) 
(219) 
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Written in explicit form, we have for i ~ I 

ui=uo + l3ihF{to,uo). (220) 

The i — 2 step yields 

W2^a2hF{to,ua) + hF{ti,ui), (221) 
U2^uo + h [(/3i + P2a2)F{tQ, uq) + p2Fiti,u,)] . (222) 

and the i = 3 step gives 

W3^h [a2a3F(io, uo) + asFih.ui) + F(i2, ^2)] (223) 

U3 = uo + [f3i + a2{f32 + 03/33)] hF{to, uq) + (^2 + f33a3)hF{ti,ui) + ^3/1-^(^2, "2) • (224) 
The corresponding times can be calculated by putting _F = 1 . This yields 

h=ta + Pih, (225) 

^2 =io + /j[/^i+/32(l + a2)], (226) 
t3 = to + h [A + /33 + (1 + a2)(/32 + 03/33)] • (227) 
The last expression can also be written in the form 

t3 = to + h {/3i + /32(1 + 02) + /33 [1 + (1 + 02)03)]} (228) 

Next we need to determine the conditions that the scheme is indeed of third order. This can be done 
by considering the differential equation 

du/dt ^ It, u(0) = uo, (229) 

for u — u{t), where uq is the initial value of u. The exact solution of eq. (22£) is uqc*. Its Taylor expansion 
for t = to + h is 

u{to + h)^ uo[l + h+^h^ + ^h^ + ...]. (230) 



The solution based on eq. (|224|) is 

U3 = uo + [/3i + a2(/32 + 03/33)] huo + {P2 + 03/33)/imi + (33hu2 (231) 



In order to compare with eq. (230) we need the explicit expressions for ui and M2, which are 

ui = {l + hpi)uo (232) 
U2 - {1 + /i [/3i + (1 + 02)/32] + /i'/3i/32} UQ. (233) 

and so we can write 

M3 = "0 + 71^ + 72/1^+73^^ (234) 

with 



71 = /3i + /33 + (1 + 02)(/32 + 03/33) 




(235) 


72 = P1P2 + /33[/32(l + 02) + /3i(l + 03)] 




(236) 


73 = /3i/32/33 




(237) 


In order for the scheme to be third order we have to require 71 = 1, 72 = 
eq. (230). Thus, we have now three equations for five unknowns. We now 
two more equations to solve for the five unknowns. 


1/2, and 73 = 
have to come 


1/6; see 
up with 
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If we assume that the intermediate ti mest eps ar e ev ahiated in equidistant time intervals, we have to 
require that the time increments in eqs. (|225|) and (226) are 1/3 and 2/3, respectively. This yields 



/3i = /32(l + «2) = 1/3, 



with two particular solutions^ 



a2 = -2/3, a3 = -l, A = 1/3, /^a = 1, /Js = 1/2, 
a2 = -1/3, as = -1, Pi = 1/3, (32 - 1/2, /Sg - 1. 



(238) 



(239) 



These are in fact the simplest 2iV-RK3 schemes that also lead to comparatively small residual errors. 

Alternatively, one can move the times closer to the end time of the time step and evaluate the right 
hand side at times — to = \h and t2 — to ~ h. This gives the particular solution 



a2 



-1/4, as = -4/3, pi = 1/2, /Ja = 2/3, Pa = 1/2. 



Again, there could be other solutions. 

Another possibility is to require that the inhomogeneous equation 

du/dt = r, u(0) = 0, 



(240) 



(241) 



is solved exactly up to some n. The exact solutions for t = h are u — ^h? for n — \ and u ^ for 



The case n — wa s alr eady considered in eqs. (225)-(227). For ?i = 1 we have F{ti,ui) — and 



F{t2,U2) = so eq. (224) gives 



"3 = Wo + {P2 + P3a3.W + /?3|/l^ 



W3 = U0 + 5[/32+/33(a3 + 2)]/l2. 

Comparing with the exact solution this yields the additional equation 

/32+/33(«3 + 2) = §. 



For n = 2 we have F{ti,ui) — g/i^ and F(t2,U2) — ■, so eq. (224) gives 



4 1.3 



U3 = + (/32 + /33a3)i/j + P^lh 
or 

U3 = uo + i[/32+/33(a3 + 4)]/i^ 
Again, comparing with the exact solution one obtains 

/32+/33(a3+4) = 3. 

This gives the solution 

a2 = -17/32, as = -32/27, /3i = 1/4, P2 = 8/9, /Jg = 3/4, 



(242) 
(243) 

(244) 

(245) 
(246) 

(247) 
(248) 



which implies that the right hand sides are evaluated at the times ti — — jh and t2 — to — In 
tables and || this scheme is referred to as "inhomogeneous" . 

Yet another idea (W. Dobler, private communication) is to obtain the additional two equations by 
requiring that the quadratic differential equation du/dt = v? with mq = 1 is solved exactly. The solution 

thank Petri Kapyla from Oulu for pointing out the second of these solutions. 
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is M = ( 1 — t ) ^, of which we only need the expansion up to /i^, so we have ^ 1 + h + . Again, we 
use eq. (224), but now with F = u'^, 



U3 = l + [Pi+ a2{P2 + asPs)] h+il32 + a3l33)hul + Pshuj. (249) 
We need and U2 only up to the term linear in h. Using eqs. ( |232| ) and ( |233| ) we have 

ul = l + 2l3ih + 0{h^), ul^l + 2[(3i + {l + a2)(32]h + 0{h^). (250) 
Inserting this in eq. ( ^49[ ) yields 

U3^1 + Sih + S2h'^ + 0{h^), (251) 

with 

Si=l3i+ a2{l32 + asPs) + {^2 + a^p^) + fi^ (252) 

and 

h = 2Pi{p2 + a^Ps) + 2/33 [Pi + (1 + a2)P2] ■ (253) 
Thus, the two additional equations are 

Pi + {1 + a2Kp2 + asPs) + p3 ^ 1, (254) 

Pi [P2 + (1 + a3)p3] + (1 + a2)P2 = 1/2. (255) 

The numerical solution is 

a2 = -0.36726297, as = -1.0278248, Pi = 0.30842796, P2 = 0.54037434, P3 = 1. (256) 

which implies that the right hand sides are evaluated at the times ti — to = 0.308/i and t2 — to = 0.650h. 
In tables 0and| this scheme is referred to as "quadratic" . 

C Derivation of the jacobian for transformation on a sphere 

Here we give the explicit derivation of eqs. ( [79|) and (|80|). We first use the transformation in the form 

x — x'^/f, y = xy/r, iix>y, (257) 

X = xy/f, y — y'^/f, H x < y. (258) 
To obtain the jacobian we differentiate with respect to x, so 

^ ^dlnx dlnf ^ dlnx dlny dlnf _ ^ _ f259) 
dhix 9 In a;' 9 In a; 9 In a; i91na;' ~ ^' 

dlnx dlnij dlnr dlnij dlnr „ , , 

1 = 77^ + ^-7^' = 2—^--—, ifx<y. (260) 
a m a; o mx omx a In a; a In a; 

We now differentiate with respect to y, 

9 In a; 9 In f 91na; 91ni/ 91nf „ , , 

= 2^^^ 1 - + - 7T^, iix>y, (261) 

91ny 91ny omy 9 my omy 

d\nx d\ny dlnf „91ny 91nf . 

0=1T1 + Ti ' 1 = TT, — , iix<y. (262) 

9 In y 9 In y omy 9 In y 9 m y 



The derivatives of f can be written as 



91nf /a:\^91na: /y\^9 In y 



9 In a; \ f J 9 In a; \f J 9 In a;' 
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(263) 



d\nf f d\nx f d\ny 



dlny \f J dlny \f J dlny 



In all cases we have 



so 



X X 

y 



^ 9 In i 9 In y 
9 In a; (9 In x ' 
9 In i 9 In y 



- 1 = 



9 In y 9 In t/ ' 



and so 



so 



^^^ ^Inx fx\^ d\nx fy\^f d\nx ^ iix>~ 
d\nx \f J d\nx \f ) \d\a.x ^ ^ — Vi 

^ dlnx dlny /i\^91n;r f y\^ f dhix ^. 

9 In a; 9 In a; \r/ 9 In a; \r/ \91na; ^ — V^ 

d\nx fy\^ .„ - 



dim dh^^_(l)\ ,is:>y, 



and so 

i91na; V^/' i91na; 

and correspondingly 

^ 9 In ai 9 In y / a: \ ^ 9 In .r / y \ ^ / 9 In a; ^. . ^ _ ^ _ 
9 In a; 9 In a; \^/ 9 In a; \f ) \91nx — 2/' 

„91na; /a;\^91ni /w\^/91na; .._ 
= 2- U 4 1 , i{x<y, 



9 In a; \fj 9 In a; \fj \91n 



\.x 



so 



and so 



91nj/ /y V -r - ^ - 



9 In 



a; 



dlnx \r J dlnx 

Hence note that there is a discontinuity of the jacobian along the diagonals. Now for the 

we have 

^ 291nx f ''\ Sinai f y\ fdlnx \ , ^ _ 

91ny \f ) dlny \f ) \dlny y' — 2^' 

^ 9 In a; dlny 91na; /91na; \ _ ^ . 

91ny 91ny \f J dlny \f J \91ny /' ^ — 2/' 

91n5 / yV ~ ^ ~ 
dlny \r J 

dlnx ^ ^ / yV dlny ^ ^ yV if J > ~ 
dlny V/' 91ny V/' ^ — V^ 



so 



and so 
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and correspondingly 



= 



9 In X 9 In y 
dhiy d In y 

dlnx 

-1 = 2 



~ X 2 

X 

r 



d\nx 
d\ny 



~\ 2 

V 



diny 



~ \ 2 
X 

f 



d\nx 



dXny 



~\ 2 

y 



dlnx 
d\ny 

dhix 



so 



and so 



dlnx 
dhiy 



= 



= -1 



d\nx 
dlny 



1-1 ^ 

r 



~\ 2 

V 



~\ 2 

y 



dhiy 
dhiy 



dhiy 
X <y 



~\ 2 

V 



, h X <y, 
hx <y, 



ii X < y. 



(280) 
(281) 
(282) 

(283) 



So, in summary 



/dlnx dlnx 

dlnx dhiy 

9 In y 9 In y 

\ dlnx dhiy 

/dlnx dlnx 

9 In a; 9 In y 

9 In y 9 In y 

V 91na: 91ny 



-1 - 
-1 + 

-1 + 
-1 - 



~x 2 

y 



~ X 2 

X 



~ \ 2 
X 



r ^ 

~x 2 

y 
f 



+1- 
+1 + 

-1 + 
+1- 



xy\ 



r 

~'\ 2 

y 



\ix>y, 



(284) 



li X < y. 



(285) 



D Derivation of the incremental jacobian for second derivatives 



Here we present the explicit derivation of eq. (110). To calculate the second derivative of a function / 
that is represented on a coordinate mesh i, is given by 



dxjdxi 



d f df 



d f df dxk 



dxi \ dxj J dxi \ dxk dx 

^ 9^ dxi dxk ^ df d^Xk 
dxidxj dxidxk dxi dxj dxk dxidxj 



or 



d^f d^f ^ J 

dxidxj dxpdxq dxk 



(286) 

(287) 
(288) 



which is just eq. (|108|), using eq. (109) for the definition of Kkij of the second order jacobian. 

To obtain the second order jacobian by successive tensor multiplication we differentiate twice the 
evolution equation for a?, 

(289) 



4"^ + nr St, 



so 



2^(n) 



d X 



2„,(n) 



d u 



dx. 



{n+l)Q^(n+l) 



9a;,- dx 



dx. 



St. 



(290) 



The expression on the left hand side is just the derivative of a Kronecker delta, see eq. (102), so it is zero. 
Thus we have 







92xi") 



9 



(") 



dx 



dx 



(ri+l) 



9ul"^ dx 



^94"^ dxf+^\ 



St, 



(291) 
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a2 (") / a (") \ r)2 (n) 



which can be written as 

= M,,Kg+zifc,p,jS^ji"^Jt, (293) 

so 

Kg=-<5i(M-i),,u,,4:)j(;\ (294) 

which is just eq. ( |112D . 

We now need to derive the equation that relates the incremental second order jacobians to the second 
order jacobian of the previous timestep. To this end we begin with the second order jacobian at time 
25t, so 

~ dx^dxf d:^\dx^dx^l ^ ' 

or 

«2^(0) o (1) o (1) o2„(l) aJO) 

^^(0^2) _ C a^fc OXp OXq ^ a Xq OXj. ^^gg^ 

Ox'^^'dx^q^ dx^p dx^p dx^'^'dx^p dx^q^ 

or 

^-^ = ^-''^^+^r''^- (297) 

For the next step we have 

k(o->3) ^ 9^4"^ _ J_ (dxl^Ox^dx^\ 

" dx^^dx^ dx^ \dx^^ dx^^ dxf J ' ^ ^ 



so 



^(0-3) ^ ,^(0^1) .(1) .(1) .(2) .(2) .(0^1)^(1) .(2) .(2) j(0) ,(1)^(2) , . 

kij kpq -"pr -"qs -"ri -"sj ' -"kl '^Irs-'ri sj ' -'kr'^rs sij' y^^uuj 

This can be written as 



which, for the general step from to n, becomes eq. (110) 



E Solution for a and r/t quenched a^-dynamo 



Here we present the explicit derivation of eq. ( |181[ ). According to mean- field theory for non-mirror 
symmetric isotropic homogeneous turbulence with no mean flow the mean magnetic field is governed by 
the equation 

= V X (aB - r;T77o J) , (301) 

where bars denote the mean fields and r/x = + ?7t is the total (microscopic plus turbulent) magnetic 
diffusion. Both a-effect and turbulent diffusion are assumed to be quenched in the same way, so 

a = =2 ' '^t = =2 • (302) 

1 + as S IBI^ 1 + ?/B B IBI^ 

In the following we assume = rye and denote 

aB/< - l/Sg - ?7b/< (303) 

We emphasize that only the turbulent magnetic diffusivity is quenched, not of course the total one. It is 
only because of the presence of microscopic diffusion that saturation is possible. 
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2 

In the simulations B is to a good approximation spatially uniform. Defining the magnetic energy as 
M{t) = j(b'^) wc have 

= 2M, (304) 

which is only a function of time. 

Consider the particular example where the large scale field varies only in the z direction eq. ( [30 ID 
becomes 

'^^ = -aB'y + r,TB'^, (305) 

iy^+aS'^+TjTB'y, (306) 

where dots and primes denote differentiation with respect to t and z, respectively. Since a < 0, the 
solution can be written in the form 

B,(y,t) 6, (^)cos(z + </)), (307) 

By{y,t) = by {t)sm{z + (308) 
where (t) and b^ (t) are positive functions of time that satisfy 

bx = \a\by - riTbx, (309) 

by = \a\bx - Tl^by. (310) 
We now choose the special initial condition, bx = by = b, so we have only one equation for the variable b. 
Note also that in the quenching factor B = b^ + by ^ 26^. Thus, we have 

db _ aoki - rjmkl ^ 



At 1 + 252/^2 

Multiplying with b yields 



■b-riktb. (311) 



dt 1 + 262/5^ 

Using the definition AI — ^b^ we have 



(312) 



J_dM^aofci-^,ofc2 
2M dt l + 2M/Mo ^ ^ ' 

where Mq — \B^. Thus, wc have 

where we have defined 77x0 = ??to + "H- We define the abbreviations A = a^ki — ij^okf for the kinematic 
growth rate of the dynamo and Ato — Vtokl for the turbulent decay rate if there were no dynamo action, 
and arrive thus at the integral 

f'' I + 2M'/Mo dAr_ 

1 - i2rikl/X) M'/Mo M' " ^'''^ 

We now also define the abbreviation 

M = 2M'/Mo (316) 

and have 



l + M dM 
JM,^,~{vkl/\)M M 



2 - 2Ai. (317) 
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which can be spht into two integrals, 



M 



dM 



Mi„ [1 - ivkf/X) M]M Jm,^, 1 - m/X) M 



M 



dM 



= 2\t. 



To solve these integrals we note that 



dx 



dx 



\ — x/xq 



dx 



{1 — x/xqi)x J Xo — X 
so eq. (318) becomes 



da; 

X 



-xo ln(xo - x), 

: Inx — ln(l — x/xq) = In 



xa - X 



In 



/ M 



X 



■In 



1 



\X/rikf — M J Tjkf \X/rik1 - M 
where tg is an integration constant. Exponentiation yields 

^ e2A(t-to) 

{X/r]kl - My+^/'^'^'i 

which is, in terms of the original variables, B /Bq = 2b^ j Bq = 2M/Mo — A4, we have 

-2\ l + 



B_ / Ji_ B_ 



„2A(t~to) 



The final field strength, 



Bfin = lim \B{t)l 

t — yoo 

is given by requiring the denominator to vanish, which yields 

A 



3 ■ 



Rewriting eq. (323) in terms of Bfin we have 



_2 ^ 1 + X/vki 

W I 

nn / 



7]kf 



1+A/jjfef 



„2A(t-to) 



(318) 

(319) 
(320) 

(321) 
(322) 



(323) 
(324) 

(325) 



(326) 



We can express to in terms of the initial field strength, Bi^i = |-B(0)|, and if the in itial field strength is 
much weaker than the final field strength, i.e. Bi^i <^ Bfin, then we can rewrite eq. (326), in the form 



B'/{l-B'/BLr'/^"^'=Bl,. 



,2At 



(327) 



Thus, for early times we have the familiar relation 

|B|»Binie^* (|B|«Bfi„), 
whereas for late times near the final field strength we have 

(1 - bVb^ J-(1+V.fe?) « (Bi„i/B,„)2 ^2At (1^1 



Ban), 



(328) 



(329) 
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Bi 



(330) 



where tgat = ^ ^ lii(i3fin/^ini) is the time it takes to reach saturation. If the Reynolds number is large 
we have A ^ ryfc^, so 



b' 



Bl 



fin 



1 _ g-2))fel(t-tsat) 



\B\ « Sfin), 



(331) 



which is identical to the result obtained from helicity conservation. 

Note that the solution (327) is governed by four parameters: Bini, -Bfin, A, and rik\. The latter is 
known from the input data to the simulation, i?ini and A can be determined from the linear growth phase 
of the dynamo (characterized by properties of the small scale dynamo!) and so Bf^^ is the only parameter 
that is determined by the nonlinearity of the dynamo and can easily be determined from the simulations. 
Once -Bfin is measured from numerical experiments we know immediately the quenching parameters 



as = f?B 



A 



Beg 

'b^ 



and since B^^^^/B^^ fa fcf/fci we have 



rjkikf 



(332) 



(333) 



which shows that q;b and rjB are proportional to the magnetic Reynolds number. 
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